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vary  T  (x)  from  1/4  to  3/4. 

Turbulent  reattachment  criteria  defined  as  (  /u  )  du  /dt. 

I  e 

Dynamic  viscosity  coefficient. 

Potential  doublet  strength. 

Kinematic  viscosity  .. /.  . 

Coordinate  in  the  computational  plane  whicti  emanates  from 
the  body  surface  grid  points  and  is  orthogonal  to  '  in 
the  computat ional  plane. 

Ratio  of  circumference  to  diameter  of  a  circle. 

Fluid  density. 

Similarity  parameter  in  Blasius'  series  of  Equation  lb; 
and  a  grid  transformation  coefficient  defincil  b'.'  I>;ualior; 

A .  7  . 

Inner  and  outer  control  volume  surfaces. 

Grid  transformation  coefficient  defined  by  Fr,;uation  A.b. 
Flement  of  the  viscous  stress  tensor. 

Element  of  the  turbulent  stress  tensor. 

Flement  of  the  total  stress  tensor,  ’  ,  +  ... 

1 .1  t  1 .] 

Transonic  small  disturbance  potential  doublet  turn  lion. 
Transonic  small  dlsturbanci-  potential  vortex  !u!’.  titn. 

Acceleration  parameter  for  pj;essure  iter.il  ion  civi”  !  . 

Equation  4o  defined  as  2  ,I“/(tt  +  ,  )  '  t  )  . 

[.b 

.n'.’i-cid  strear,  function. 
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Symbol 

^  Magnitude  of  the  local  flow  field  vorticity,  curl  y^;  and 

a  local  acceleration  parameter. 

u'  Acceleration  parameter  for  the  field  pressure  finite 

^  difference  equations. 

..  ,  Acceleration  parameter  for  the  surface  pressure  iteration, 

pb 

Acceleration  parameter  for  the  u  and  v  velocity  component 
uv  ^ •  j  . 

finite  difference  equations. 


Subscript 

b 

e 


f 

i 


j 


k 


Denotes  a  location  on  the  body  surface. 

Denotes  a  location  at  the  edge  of  the  boundary  layer  or 
wake  . 

Denotes  a  location  on  the  computational  far  field  boundary. 

Denotes  the  f  coordinate  location  in  the  computational 
plane  for  finite  difference  equations;  and  an  indicial 
notation  index  in  differential  equations. 

Denotes  the  r,  coordinate  location  in  the  computational 
plane  for  finite  difference  equations;  and  an  indicial 
notation  index  in  differential  equations. 

A  component  index  in  indicia!  notation. 


max , min 


Denote  the  maximum  and  minimum  values  of  a  variable. 


S  Denotes  a  value  at  the  laminar  separation  point  which 

defines  the  beginning  of  the  bubble. 


w 

Denote; 

■;  a  value  in  tlic 

wake* 

region . 

x,xx 

Denote 

dif  f  erent iat ion 

with 

respec  t 

to  X 

y 

Denote 

d i [ f  erent iat ion 

with 

respec  t 

to  y 

T";  ,r’V 

Denote 

differentiation 

with 

respect 

t  0  ’ 

r  r  f 

'  ^  ' 

Denote 

di f  f  erent iat ion 

with 

respic  t 

to  ■ 

Denott's  a  vector  quantity. 


Superset  ip t 

*  Denotes  a  quantity  in  the  transformed  pl.ii-ii'. 

I'l’iioti's  a  turbulence  liric  aver.i.'td  qviantitv. 
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cript 

Denotes  a  turbulent  fluctuating  quantity. 

Denotes  a  temporary  nondimensional  variable  definition 

Iteration  number  in  the  successive-over-relaxation 
iteration  procedure. 

Pertaining  to  n  contours. 

Pertaining  to  C  contours. 

Derivative  operator. 

Substantial  derivative  operator  (D-/Dt  =  c<-/H  +  y'V-) 
Increment . 

Summation . 

Del-operator  defined  as  i  ^  +  jl~- 
Partial  derivative  operator. 
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ABSTRACT 

Numerical  solutions  are  obtained  for  two-dimensional  incompressible 
turbulent  viscous  flow  over  airfoils  of  arbitrary  geometry.  An  algebraic 
eddy  viscosity  turbulence  model  based  on  Prandtl's  mixing  length  theory 
is  modified  for  separated  adverse  pressure  gradient  flows.  Finite  dif¬ 
ference  methods  for  solving  the  inviscid  stream  function  equation  and  the 
incompressible  laminar  Navier-Stokes  equations  are  used.  A  finite  difference 
method  for  solving  the  Reynolds  averaged  incompressible  turbulent  two- 
dimensional  N’avier-Stokes  equations  is  employed. 

The  inviscid  stream  function  equation  and  the  N’avier-Stokes  equations 
are  transformed  using  a  curvilinear  transformation.  A  body-fitted  coor¬ 
dinate  system  with  a  constant  coordinate  line  defining  the  airfoil  section 
surface  is  transformed  to  a  rectangular  coordinate  system  in  the  trans¬ 
formed  or  computational  plane.  An  elliptic  partial  differential  Poisson 
equation  for  each  coordinate  is  used  to  generate  the  coordinate  system  in 
the  physical  plane  for  arbitrary  airfoils. 

The  two-dimensional  time  dependent  Reynolds  averaged  incompressible 
Navier-Stokes  equations  in  the  primitive  variables  of  velocity  and 
pressure  and  a  Poisson  pressure  equation  are  numerically  solved.  Turbu¬ 
lence  is  modelled  with  an  adverse  pressure  gradient  eddy  viscosity  tecli- 
nique.  An  implicit  finite  difference  method  is  used  to  solve  the  set  of 
transformed  partial  differential  equations.  The  system  of  linearized 
simultaneous  difference  equations,  at  each  time  step,  is  solved  using 
successive-over-relaxation  iteration.  Far  field  boundary  conditions  are 


examined.  Solutions  for  a  NACA  0012  airfoil  at  angles  of  attack  varying 
from  five  to  11.5  degrees  at  a  chord  Reynolds  number  of  170,000  are 
obtained.  Velocity  profiles  near  the  airfoil  surface  and  surface  pres¬ 
sure  distributions  are  presented  and  compared  with  experimental  data. 

Lift  and  drag  coefficients  agree  well  with  experimental  values.  The 
computed  lift  coefficients  near  stall  are  within  57,  of  the  experimental 
measurements ,  and  the  numerical  drag  coefficients  agree  within  ten  drag 
counts  in  the  region  of  maximum  lift  to  drag  ratio.  The  short  laminar 
separation  bubble  near  the  suction  pressure  peak  is  numerically  determined. 
The  variation  of  bubble  length  and  turbulent  transition  length  with  angle  of 
attack  are  similar  to  experimental  trends. 
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SECTION  I 


INTROnrCTION 


Much  effort  has  been  expended  by  the  aeronautical  coriruin.itv  in  di  - 
ternininp  the  aerodynamic  characteristics  of  airfoils.  1. inear  c.etisu:. 
are  extensively  used  in  design  work  for  studying  configurations  at  sr..i]l 
angles  of  attack  with  negligible  flow  separation.  Experimental  wind 
tunnel  investigations  are  used  to  determine  the  charactcrist ics  near  stall 
where  separation  phenomena  become  important.  Recent  developments  in 
nvimerical  techniques  have  stimulated  research  on  another  approach,  namelv 
the  nuraerical  solution  of  the  Kavier-Stokes  equations.  These  equal  ion.s 
model  the  vi.scous  effects  which  contribute  to  airfoil  stall.  For  tl'.is 
reason,  numerical  Navier-Stokes  methods  offer  the  possibility  of  determin¬ 
ing  the  aerodynamic  characteristics  for  airfoil.s  experiencing  stall. 
Numerical  methods  can  also  complement  experimental  methods  by  efficiently 
extending  the  range  of  parameters  under  investigation.  Furthermore, 
numerical  methods  eliminate  model  support  interference  and  wall  interference 
effects  found  in  wind  tunnel  testing. 

The  purpose  of  this  investigation  is  to  develop  a  numerical  Navier- 
Stokes  method  that  will  accurately  determine  the  aerodynamic  chia  rac  t  er  i  s',  i  es 
of  incom.pressiblc  turbulent  viscous  flow  over  two-dimensional  airfoils  lu'.ir 
stall. 

The  development  of  a  nuraerical  method  for  turbulent  flow  requires 
a  suitable  technique  for  distributing  points  throughout  the  flow  !  i  i  1  d. 
and  a  model  which  describes  the  behavior  of  turbulence  within  sjnci'ie 
regions  of  the  flow  field.  A  survey  of  numeric.il  grid  general  in,-  tee'nj; 
and  available  turbulence  models  i.s  presented.  The  (juanlit'.  of  1  i  1 1  r.i :  i;  ta 
conct'rncv!  with  the  N.iv  i  er-St  okes  equations  is  extinsive.  I'bi- ta  '  !>  .  .i 
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summary  of  the  literature,  which  describes  formulations  of  the  Navier- 
Stokes  equations  and  their  numerical  solvint;  techniques  applied  to  flows 
over  airfoils,  is  given.  The  research  objective.s  for  tliis  work  ari- 
then  discussed,  and  a  summary  of  the  remaining  sections  i.s  presented. 

Body-fitted  curvilinear  coordinate  systems  greatly  enliance  tl\i' 
application  of  numerical  methods  to  practical  boundary  value  proble::.s 
involving  partial  differential  equations.  The  representation  of  a 
boundary  surface  as  a  coordinate  line  reduces  the  difficulties  assoc i.ated 
with  numerically  specifying  boundary  conditions  by  interpolation  ii'.  finite 
differences  methods.  In  the  physical  (x,y')  plane,  values  for  oiu'  cov.’.pc.- 
tational  coordinate  f  are  .specified  at  selected  points  on  hoth  tiie  ho.lv 
surface  and  the  outer  boundary.  Constant  value.s  for  the  other  com.’pvitat  iona 
coordinate  r  are  specified  on  both  the  body  surface  atid  tlie  outer  flow 
boundary.  The  transformed  computational  ’)  plane  tlicn  becomes  .1 
rectangular  region  with  an  orthogonal  grid.  Winslow  (1)  and  Chu  (.2) 
introduced  the  concept  for  two-dimensional  regions  interior  to  a  clos<.i 
boundary.  Their  transformed  coordinates  are  solutions  of  Laplace's 
equation  in  the  physical  plane  and  define  a  triangular  mesh  systeti  in 
the  physical  plane.  Amsden  and  Hirt  (31  took  the  physical  coordinates 
to  he  .solutions  of  a  modified  Laplace's  equatioti  in  the  t  ransf  orr.cd 
plane . 

Thompson,  Thames,  and  Mastin  (s,b,bl  generalised  tlic  method  for  iLe 
auton.itic  generation  of  body-fitted  coordinates  for  any  two-d  imeiwv  i  on.!  1  , 
mu  1 1  i -connec  ted  region.  They  also  introdnci'd  the  use  of  forcing  fi:i:ct  i.'n.- 
in  a  Poisson  equation  wiiich  proviiies  nesh  contrti]  in  reg.ioiis  with  larg.e 
gra.iients.  Hodc.e  (71  developed  an  .nuomated  grid  line  attraction  r.ith.'d 
based  en  houndarv  laver  theory  which  determinc-s  the  coe  f  f  i  c  i  ..ni  t  s  i-'  iii,-. 


forcirii;  function.  Hodge  assuned  a  Blasius  houndary  lavL-r  profiii.  arid 
distributed  his  grid  points  at  approximately  equal  velocity  i n c r eraT. t ^ 
in  the  boundary  layer.  Steyer  and  Sorenson  (S)  introduced  auxiliary 
conditions  for  the  forcing  functions  which  provide  ancle  and  distance 
control  at  the  inner  boundary  surface.  liie  angle  with  w'nich  c.  1  i;u 

intersects  the  body  surface  is  specified  by  a  functio”.  f  )  ,  an.d  the 
rate  of  change  of  arclength  with  ’  on  a  f  line  at  th.-  bodv  surface  i 
prescribed  by  s,  (  0  .  Sorenson  (9)  later  Imposed  sir.ila.r  condition.c  on, 
an  outer  computational  boundary.  Tl'.e.se  auxiliary  conditions  are  used  to 
solve  for  coefficients  in  the  cho.sen  exponential  forcing  functions.  he 
geometric  conditions  hold  exactly  only  in  the  limit  as  ’  approaches  ccro. 
Sorenson  reported  that  numerical  instahil  it  ies  occur  for  large  ch.a.n.cc. 
in  the  coef  f  icitints  during  successive  iterations  and.  for  boimdaries  v.ith 
sharp  corners.  He  implemented  a  lim.it  function  which  damped  the  ch.inge  of 
the  value  for  each  coefficient  over  one  iteration;  and  at  sharp  corners  ik 
computed  average  values  of  each,  coefficient  from  data  at  the  neighbor iiy 
boundary  points,  ’'.astin  and  Thomp.son  (101  have  also  extendid  th..  1  1  i  I'f 
body-fitted  coordinate  generation  technique  to  three  dimen.sions  for  .sir:y  lL 
geomet  r ies  . 

Some  other  specific  grid  generation  techniques  wb.ich  usv  properties  of 
elliptic  differential  equations  have  also  been  reported.  Meyde r  (111 
constructed  an  orthogonal  curvilinear  coordinate  systec:  by  usir.g  Llectric 
field  theory.  He  solved  the  potential  equation  twice  witii  different  sets 
of  mixed  Dirichlet  and  Neumann  boundary  conditions  for  the  electric  ]iotenli,!l 
and  electric  f^'rce  lines,  respectively.  Th.ese  soUitions.  h.owLver,  wi  i'l' 
obtained  in  the  physic.il  plane  using  in  t  erpola  1 1  on  and  becaru  iIk  cur-- i  1  i-.  .  v 
coordinate  line.s  in  the  physical  plane.  Tin  coordin.'.te  r-.etrics  w..  r.  th.:, 
used  to  formul.-itt-  a  finite  difference  equation  ■w}:;cli  w.js  als.  .^cl\-id  in  Ih. 
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phvsical  plant'.  Conformal  mapping  techniques  vliich  employ  the 


1 heoJor^en- 

Garick  (11)  trans:  -rmation  have  been  examined  by  Ives  (13i.  he  introduced 
the  use  of  fast  Fourier  transform  methods  and  developed  a  new  class  of 
transf  ormat  ions  which  naps  the  flow  field  of  a  two-eler.ent  airfoil  or.to 
the  regioj;  between  two  concentric  circles.  Conformal  techniques  are  not 
capable  of  being  extended  to  three-dimensional  geometries  (-).  heither 
of  tfiese  approaches  offers  a  convenient  means  for  mesh  control  in  regions 
of  large  flow  gradients. 

Kc'Centl;,-,  a  geometric  grid  generation  technique  has  been  introduced  by 
Gibe-ling,  Sbiam.roth,  and  Eiseman  (14).  The  technique  param.eter izes'  (t)  the 
body  and  outer  boundary  surfaces  and  uses  a  stretching  function  H(r}  for 
mesh  attraction  along  constructed  lines  perpendicular  to  t!ie  body  surface. 
I'nit  increments  for  ordered  pairs  (t,r)  generate  the  corresponding  computa¬ 
tional  plane.  Further  refinements  which  provide  angle  and  arclengtb. 
variation  control  at  a  body  surface  have  subsequently  been  developed  b;c 
Eiseman  (13,16,17)  . 

The  searcii  for  an  accurate  and  universal  turbulence  model  has  paralltled 
the  development  of  body-fitted  grid  gener-tior;  techniques.  In  principle, 
the  N'avier-Stokes  equations  completely  describe  the  turbulent  f  luc  tuat  i:-.g 
fluid  motion.  The  required  mesh  resolution,  necessary  to  resolve  the  tur¬ 
bulent  eddies  with  varying  length  scales,  when  translated  into  cc'::.:  :ttr 
resources  presently  make  this  approach  unfeasible.  Many  quantitio  of  im- 
gineering  Interest  in  turbulent  flows  involve  a  mean  value  take-r.  over  a  time 
interval.  The  time  interval  is  sufficiently  long  to  include  man"  f 1 uc I uu t  ions 
while  small  compared  with  the  characteristic  time  of  tbje  me<ni  Ih",;.  lie 
Mavier-Stokes  equations  can  then  be  re- f ormulated  using  these  r.e.;n  flov. 
variables.  This  I’.eynolds  averaging  procedure  introduees  the  -a  a",  of  te  rm-; 


involving  fluctuating  quantities. 


The  Reynolds  stress  components 


-V  u!u!  are  the  most  common  terms  of  this  type.  In  order  to  solve  the 
1  J 

a.’eraged  form  of  the  Navier-Stokes  equations,  "turbulence  closure"  must  be 
achieved  by  suitably  modelling  these  additional  terms.  This  approach  to 
turbulence  has  led  to  models  which  introduce  auxiliary  relationships  ranging 
from  algebraic  equations  to  several  partial  differential  equations.  Tiiese 
models  are  commonly  categorized  and  are  now  sumjnarized  according  to  trie 
number  of  additional  partial  differential  equations  which  comprise  the  miodel. 

The  algebraic  or  zero  equation  models  have  their  origins  in  Boussinesc's 
(18)  eddy  viscosity  hypothesis  and  Prandtl's  (19)  m;ixing  length  model.  The 
local  turbulent  stresses  are  assumed  proportional  to  the  local  mean  flow 
strain  rates  with  the  proportionality  constant  defined  as  an  equivalent 
or  eddy  viscosity.  The  eddy  viscosity  m.odels  of  Ceheci-Smith  (2Q)  ,  I'.ellor- 
Herring  (21),  and  Patankar-Spalding  (22)  represent  this  approach.  The 
boundary  layer  in  each  method  is  divided  into  an  inner  near  wall  layer  and 
an  outer  wake  layer  with  separate  expressions  for  the  eddy  viscosity  co¬ 
efficient.  In  the  Cebeci-Smith  and  Patanker-Spalding  models  the  inner  mixing 
length  varies  as  a  linear  function  of  the  normal  distance  from  the  wall  modi¬ 
fied  with  a  Van  Driest  (23)  laminar  sublayer  correction.  The  Cebeci  and  Mellor 
models  both  base  the  outer  length  scale  on  the  displacement  thickness  while 
the  Patanker-Spalding  model  uses  the  boundary  layer  thickness  directly. 

These  models  have  been  successfully  extended  and  applied  (24)  to  a  wide  variety 
of  boundary  layer  flow  geometries  involving  compressibility,  heat  and  mass 
transfer,  and  curvature  effects.  In  addition,  the  mean  flow  models  are 


computationally  efficient.  Launder  and  Spalding  (25)  point  out,  however, 


that  these  models  predict  a  vanishing  eddy  viscosity  where  the  velocity 
gradient  is  zero  and  have  not  been  successful  for  large  separated  recirculat¬ 
ing  flows. 
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The  one  equation  turbulence  model  introduced  by  Prandtl  (26)  is  an 


extension  of  the  algebraic  technique.  In  this  approach  the  solution  of 

the  partial  differential  turbulent  kinetic  energy  equation  usually  provides 

a  local  velocity  scale  given  by  q'  =  ulu!^  where  suinmation  is  implied. 

The  length  scale  L  is  prescribed  by  an  algebraic  expression  as  before.  Terms 

involving  fluctuating  quantities  other  than  q  must  still  be  modelled. 

Glushko  (27),  Mellor  and  Herring  (28),  and  Wolfshtein  (29)  model  the  gradient 

diffusion  term  -{(u'^  p/.  )  +  u'j^  ^'^'i  '^'i  with  a  gradient  of  q  given 

as  (q'/2)  where  is  a  specified  constant  or  function  and  ■  is 

Q  k  Q 

the  eddy  viscosity.  The  Reynolds  stress  is  related  to  the  mean  flow  as 
in  the  algebraic  approach  except  that  the  eddy  viscosity  is  assumed  propor¬ 
tional  to  qL.  Bradshaw,  et  al  (30)  model  the  gradient  diffusion  tern,  with 
an  expression  G  q”  where  G  is  an  empirical  constant  or  function  and 
is  a  velocity  characteristic  of  large  eddy  motions.  They  also  write  the 
Reynolds  stress  directly  as  a  function  of  q  in  the  xorm  T  ,,=A,.  q'.  In 

tij  ij  ^ 

each  model  the  isotropic  dissipation  is  modelled  by  a  form  q  7L.  Nee  and 
Kovasznay  (31)  use  a  rate  equation  for  the  total  viscosity  N'  =  i 
in  place  of  the  kinetic  energy  equation.  They  formulated  expressions  for 
the  generation  and  dissipation  terns  involving  various  constants  and 
the  length  scale.  Launder  and  Spalding  (25)  observe  that  the  one  differential 
equation  models  require  a  moderate  increase  in  computer  resources  but  do 
not  in  general  provide  more  accurate  results  than  the  results  obtained  from 
algebraic  methods. 

Kolmogorov  (32)  introduced  the  somewhat  more  complex  type  of  turbulence 
model  which  uses  two  partial  differential  equations.  A  form  of  the  turbulent 
kinetic  energy  equation  provides  a  local  velocity  scale.  The  local  length 
scale  is  obtained  from  a  second  equation.  Kg  and  Spalding  (33)  formalized 


6 


this  approach  by  introducing  the  energy-length  equation  derived  by 
Rotta  (34).  They  also  used  the  Glushko  closure  model  in  the  turbulent 
energy  equation.  This  '.lass  of  two  equation  models  is  named  the  k  -  kL 
turbulence  model  where  the  eddy  viscosity  is  proportional  to  k‘  L.  Saffman 
(35)  has  used  a  transport  equation  for  the  mean  vorticity  together  with 
the  turbulent  kinetic  energy.  The  eddy  viscosity  becomes  proportional  to 
k/-^  in  this  k  -  -  class  of  closure  models  for  turbulence.  The  local  length 

J 

scale  is  assumed  to  be  proportional  to  k •  Wilcox  and  Rubesin  (36)  h.avt- 
modified  this  approach  for  compressible  flows  and  generalized  the  constit¬ 
utive  equation.  Jones  and  Launder  (37)  use  a  transport  equation  for  tlio  rat 
of  dissipation  of  turbulent  energy  g  along  with  the  turbulent  kinetic 
energy  equation.  Glushko  type  closure  models  are  assumed  for  the  terms  in¬ 
volving  fluctuating  quantities  in  both  equations.  This  k  -  -t  class  of  two 

3/2 

equation  models  has  a  local  length  scale  proportional  to  k  "/'  with  the 
eddy  viscosity  proportional  to  k  /-:  .  The  two  equation  models  have  been 
applied  to  various  boundary  layer  and  free  shear  layer  flows  with  a  variety 
of  constants  and  closure  models.  The  two  equation  models  require  signifi¬ 
cant  increases  in  computational  resources  and  have  not  led  to  a  model  of 
universal  applicability  (24,25). 

The  search  for  a  more  general  turbulence  model  has  led  to  the  use  of 
the  transport  equations  for  the  Reynolds  stresses.  In  this  approach  the  edd 
viscosity  concept  is  discarded.  Closure  models,  however,  are  still  requirid 
for  the  terms  containing  fluctuating  quantities  other  than  the  Reynolds 
stresses.  Donaldson  (38)  introdu.  closure  models  which  express  the 
fluctuating  terms  as  functions  of  the  Reynolds  stress  and  chosen  lenglli 
scales.  Hanjalic  and  Launder  (39)  introduce  closure  models  which  use 
the  Reynolds  stress  and  retain  the  equations  for  turbulent  kinetic  energy 
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Only  cases  with  the  siinplified  two- 


k  and  turbulent  dissipation  . 
dernensional  boundary  layer  approximations  have  been  investigated.  The 
extensive  computational  resources  and  initial  state  of  development  of  these 
models  preclude  this  type  of  approach  from  practical  consideration  as  a 
turbulence  model  for  a  complex  flow'  field  computation. 

Theoretical  turlilence  models  of  further  complexity  have  appeared. 

Kolovandin  and  Votutin  (40)  introduced  a  statistical  theory  where  additional 
equations  are  obtained  for  other  correlations  of  the  fluctuating  quantities. 
Ferziger  (41)  took  a  meteorological  viewpoint  of  turbulent  flows  by  numeri¬ 
cally  simulating  large  scale  eddies  while  modeling  the  small  scale  structures 
with  an  eddy  viscosity  technique.  This  approach  is  a  first  step  toward 
numerically  solving  a  turbulent  flow  field  with  the  instantaneous  equations 
of  notion. 

The  review  of  available  turbulence  models  provides  the  basis  for  select¬ 
ing  a  suitable  approach  for  use  in  the  numerical  solution.  The  algebraic 
and  the  two  equation  eddy  viscosity  iDodels  appear  to  be  realistic  choices. 

As  previously  reported  (25),  the  one  equation  techniques  yield  unimproved 
results  compared  with  algebraic  methods  and  at  additional  cost  in  computer 
resources.  The  two  equation  methods  require  the  solution  of  additional  partial 
differential  equations.  In  these  methods,  terms  involving  fluctuating  quanti¬ 
ties,  except  the  Reynolds  stresses,  must  still  be  modelled  using  additional 
coefficients.  For  these  reasons,  the  simpler  algebraic  eddy  viscosity  technique, 
which  requires  considerably  less  computer  resources,  is  employed.  If  the 
physical  phenomena  associated  with  separated  adverse  pressure  gradient  flows 
can  be  included  in  the  turbulence  model,  then  the  accurate  calculation  of  the 
aerodynamic  characteristics  may  be  accomplished  with  an  algebraic  technique. 

A  survey  of  the  previous  research  involving  numerical  solutions  of 
the  N'avier-Stokes  equations  for  flow  over  airfoils  establishes  the  prc- 
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diction  level  of  computational  methods  and  also  reviews  the  aumerical 
algorithms.  The  grid  generation  technique  and  turbulence  model  employed, 
when  applicable,  are  also  included. 

Early  numerical  solutions  of  the  Navier-Stokes  equations  for  flow  over 
airfoils  jsed  the  vorticity-stream  function  formulation  with  an  automated 
grid  generation  technique  (4).  Walker  (42)  applied  the  metViod  to  the  laminar 
flow  over  a  flat  plate  and  compared  the  numerical  solution  with  Blasius' 

(43)  solution.  Thames  (44)  used  body-fitted  coordinates  with  the  vorticity- 
stream  function  approach  and  solved  the  Navier-Stokes  equations  for  various 
bodies  in  doubly  connected  regions.  He  obtained  solutions  for  the  flow  over 
airfoils  at  chord  Reynolds  numbers  less  than  10  .  Problems  mainly  attributed 
to  wall  vorticity  developed  for  solutions  of  airfoils  at  angle  of  attack. 

>!ehta  and  Lavan  (45)  also  used  a  vo'  ticity-stream.  function  meth.od 
in  studying  the  laminar  starting  vortex  and  separation  bubbles  for  impul¬ 
sively  started  incompressible  laminar  flow  over  a  Joukowski  airfoil  at 

A 

Reynolds  numbers  less  than  10  .  They  used  three  point  backv.’ard  time 
differences  and  centered  spatial  differences.  The  computational  grid 
was  obtained  through  a  conformal  transformation  followed  by  a  radial 
stretching  transformation. 

Reddy  and  Thompson  (46)  applied  an  integro-differential ,  vorticity- 
velocity  field  method  for  the  solution  of  incompressible  flow  in  doubly 
connected  regions.  Backward  time,  centered  spatial  (BTCS)  differences 
were  applied  to  the  N’avier-Stokes  equations.  The  difference  equations 
were  solved  using  successive-over-relaxation  (SOR)  iteration.  I’u-y 
also  employed  the  coordinate  mesh  attraction  technique  with  a  time  depend.,  r.l 
expanding  mesh  system.  Symmetric  airfoils  at  zero  angle  of  attack 
with  a  Reynolds  number  less  than  10^  were  considered.  At  the  liiglur 
Reynolds  numbers,  the  calculation  of  surface  vorticitv  required  a  lar.u  nun- 
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her  of  grid  points  on  the  surface  whitli  greatly  increased  the  r<.':>;  i  rtd 
computer  time.  A  steady  state  solution  was  not  obtained. 

The  vorticity-velocity  field  formulation  has  been  appliid  by  Sanhar 
and  Wu  (47)  to  the  ca.'-e  of  incompressible  laminar  flow  about  an  o.-,ei  1  lat  in. 
airfoil.  They  used  a  12’  tluck  Joukovski  airfoil  at  a  Keyneld.s  nu:  hi  r  of 
'  1000.  Triangular  finite  elements  were  constructed  near  the  airfoil  surface 
and  a  rectangular  mesh  was  used  awav  iron  the  airfoil.  A  conform;il  trans- 
fori.iation  was  used  to  map  the  airfoil  in  the  physical  plant-  OTit.'  a  unit  circle 
in  the  computational  plane.  Sankar  and  Tassa  (48)  investigated  this  sa:'i. 
problem,  wit'u  a  compressible  flow  formulation  of  the  bavier-Siokes  e'.u.i  t  i  en  -  . 
Thiey  u.secl  a  conformal  t  ran.sf  tirm.it  ii>n  followed  by  an  algebrait  radial  stretchir.. 
transf  or::,at  ion  .  The  a  1 1  erna  t  ing-d  irec  t  ion-im;p  1  ic  i  t  (ATT)  finite  difference 
method  of  Briley  and  Mcl'onald  (49)  was  used  to  obtain  solutit'r.s  for  ••eyr.old.^ 
numbers  less  than  10  and  a  Mach  number  of  h.2.  In  a  separate  research 
effort,  Sugavanati  and  Wu  (50)  attempted  to  use  a  two  equation  k-  turbu¬ 
lence  model  with  a  vortic  ity-vt  locity  formulation.  .-5  conformal  transform.a- 
tion  for  a  12  percent  thick  .loukowski  airfi.'il  was  ust-d..  They  ex.per  ie::Ced 
difficulties  in  obtaining  a  converged  sc'lutioi;  even  ."'fier  eight  hours  of 
CYBMr.  74  CPr  time  were  expended.  Variation,-,  for  botli  lift  and  drag  of  the 
order  of  50  percent  occurred. 

The  use  of  the  primitive  variable.s  of  vcliHity  and  pres.sure  for 
incompressible  flow  was  introduced  by  Harlow  and  Wklc),  (hi)  in  tin. 
explicit  forward  time,  centered  spatial  (FTfs)  M.irk-.ind-(.  e  1 1  r,eth..l. 

Tliey  included  a  Poisson  equation  for  pressure  whi(l.  is  obtained  ly.  taking 
the  divergence  of  the  momentum  equation.  Ihc;.  .i  1  s  .  foiii  d  that  ths-  velocity 
divergence  term.s  in  the  momentum  equation^  w<  ri  r(.-;uired  for  fhi  pn-s-iirt 
field  calculation.  Hirt  and  Harlow  (hfi  iurther  developed  t  hi  r.cth'd. 

Hodge  (53)  considered  the  case-  of  lam.in.ir  ini  oc.p  re-'-,  s  i  M ,  viscon  t  1  o-..  i.t 


an  airfoil  at  angle  of  attach.  Hi-  used  .i  forr  -  f  t'.i  I' 


'ii  ,  it  .1  1  (  ) 
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body-fitted  grid  transformation  and  applied  an  implicit  BTCS  d  i  f  f  erenc  itif' 
method  to  the  Navier-Stokes  equations.  Tiie  system  of  difference  ecjuat  ions 
was  solved  with  SOR  iteration.  Various  method.s  of  calculating  tht  pressure 
field  were  investigated.  Hodge  concluded  that  the  momentum  equations  should 
retain  the  velocity  divergence  terms  and  that  a  Pols.son  pressure  equation 
should  be  used  to  satisfy  continuity.  Ghia,  Hankey,  and  Hodge  (34)  applied 
a  Poisson  pressure  equation  and  the  primitive  variables  form  of  t!ie  l.'avier- 
Stokes  equations  to  study  incompressible  driven  flow  in  a  square  cavity 
for  Reynolds  numbers  under  1000.  They  used  an  al ternat ing-direc t ion-in^pl i c i t 
(ADI)  finite  difference  technique  for  the  momentum  equations  and  SOR  iter.it 
for  the  pressure  equation.  A  Neumann  boundary  condition  derived  from  th'- 
normal  com.ponent  of  the  momentum  equations  was  employed  to  compute  the  wall 
pressure . 

The  Poisson  pressure  equation  has  been  further  examined  by  Chien 
(55)  for  internal  flows.  He  observed  that  the  calculation  of  pressure 
by  direct  integration  of  the  momentum  equations  can  be  inaccurate  when 
large  velocity  variations  are  present.  He  found  that  forms  of  the 
pressure  equation  were  more  suitable  for  computing  the  pressure  field 
near  boundaries.  Recently,  Hodge,  et  al  (7)  applied  an  implicit  backward 
time  finite  difference  method  with  both  upwind  and  centered  spatial 
differences  to  the  Navier-Stokes  equations  in  primitive  variables  form. 

A  Poisson  pressure  equation  was  used  to  obtain  a  solution  for  laminar 
flow  over  airfoils  at  angle  of  attack  with  a  Reynolds  number  of  order  If'^. 

The  solution  contains  a  large  oscillating  separated  region  whicli  gives 
approximate  trends  in  lift  but  very  poor  agreement  wit!)  a\'ail.ible  drag  data 
(5G)  . 

Parly  nuratrical  i;avi'  r-Stokes  investigations  of  turbulent  flow  c'Ver 
airfoils  considered  two-dimensional  transonic  viscous  flow  at  small  angK'. 
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of  attack.  Deiwert  (57)  applied  MacCormack  s  explicit  method  (5fe)  with 


a  50x38  rectangular  based  exponentially  stretched  mesh.  He  used  20  uni¬ 
formly  spaced  points  to  define  the  upper  airfoil  surface  and  imposed  the 
Neumann  boundary  layer  condition  hp/Nn  =  0  on  the  surface.  Steger  (59) 
used  the  implicit  Beam-Warming  method  (60)  with  a  71x33  grid  obtained  fror. 
a  modified  Thompson  transformation  (4).  He  calculated  the  solution  on  t!ic 
branch  cut  with  a  linear  extrapolation  procedure  and  evaluated  surface  pres¬ 
sures  using  the  momentum  equations  for  the  direction  normal  to  the  surface. 
Walitt,  et  al  (61)  applied  the  first  order  method  of  Irulio  (62)  coupled  wit': 
a  boundary  layer  technique  used  near  tlie  airfoil  leading  edge.  They  usid  .i 
130x68  mesh  system  but  did  not  capture  the  full  suction  pressure  on  the  upper 
surface  near  the  airfoil  nose.  Each  investigation  used  a  Cebeci-Smith  ('3) 
type  algebraic  eddy  viscosity  model  and  the  Reynolds  averaged  compress il' K 
Navier-Stokes  equations. 

Recently,  Shamroth  and  Gibeling  (64)  used  the  Br  i  ley-V.chona  1  d 
(49))  Implicit  finite  difference  formulation  with  a  constructive  tvpe 
81x30  grid  system  to  compute  turbulent  flow  over  an  airfoil  at  an  an.-li 
of  attack  of  6  degrees  and  Reynolds  number  of  10^.  They  first  tried 
a  tw’o  equation  k-e  turbulence  model  but  experienced  convergence  prol  lems 
near  the  leading  edge  and  far  field  flow  regions.  They  reported  obtaining 
large  or  negative  values  for  the  turbulent  viscosity  in  essent  i)il  I-.' 
laminar  regions.  The  turbulent  kinetic  energy  equation  with  an  a  1 ee!  ra i ea 1 1  ■' 
prescribed  length  scale  was  then  used.  Major  discrepancies  occurred  for  tin. 
mean  pressure  distribution  solution  in  the  suction  peak  and  trailing  edgv 
regions  of  the  airfoil  surface.  They  also  noted  that  the  construct itw 
grid  technique  caused  a  crossing  of  grid  lines  of  the  same  family  f(>r 
highly  cambered  airfoils. 

The  examination  of  existing  numerical  solutions  of  the  Navier-btL'l  is 
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equations  for  two-dimensional  flow  over  airfoils  reveals  Reynolds  nur.bc-r 
and  angle  of  attack  restrictions.  As  a  result,  numerical  Kavier-Stokes 
methods  for  the  accurate  computation  of  the  flow  field  and  resulting 
aerodynamic  characteristics  for  an  airfoil  near  stall  are  unavailable.  The 
development  of  such  a  numerical  method  for  solving  incompressible  tvo-dimen 
sional  turbulent  flow  over  airfoil  sections  near  stall  does,  therefore, 
constitute  a  significant  contribution  in  computational  aerodynamics.  The 
fulfillment  of  this  research  goal  requires  tiie  formiilation  of  an  adequaLe 
turbulence  model  for  use  in  botli  leading  and  trailing  edge  separated  region 
The  far  field  boundary  conditions  for  incompressible  flows  must  also  !e 
examined . 

The  selection  of  a  suitable  numerical  approach  for  incomrress ill e 
turbulent  flows  is  based  on  sever.il  considerations.  An  implicit  tec  hn  i 
is  preferred  because  of  the  required  small  grid  .spacing  nt'ar  the  body  sur¬ 
face  necessary  to  resolve  viscous  stresses.  Stability  criteria  of  explicit 
method.s  would  impose  a  small  time  step  restriction  as  a  nsult  oi  the  .;rid 
size.  Ihe  nur.erical  ru'thod  should  be  capable  of  predicting  laminar  flow 
for  Reynolds  nurr.bers  approaching  turbulent  fli'w  conditions  sin.ce  regions  of 
laminar  f  1  ow  r.,iy  occur.  The  usi  of  primitive  v.iriahli-s  is  re.adil;'  ex  t  oii.led 
to  threi  dimensuais  and  provides  for  tin  direct  L  iKul.ition  o‘  the  pris..;r. 
field.  These  criteria  are  satisfied  by  an  implicit  finite  difKrcinc  prc- 
ced'irc  developed  b’.  hedge  7)  which  uses  primitive  variai  Ic  and  sisec 

over  -  re- 1 .)  xa  t  i  on  iteration.  The  in.plicit  Mrite  difference  i: ,  ;  n  i-  i  ;  In 
in  thi-^  invi-sl  i  g.'t  i  "n  . 

TliC  rernaininv  s>cti(--ns  of  t  h  i '•  reseaieh  eftor;  dis,  ,  t  tia  tc  n;  i'n. 

T.p' 1  o-,-id  in  ach  i  e'.' i  no  the  def  in<'d  rese.ircb  .’0,1;  .  Xi  ,  ;  ii  ■,  ;  ]  ,,  ;  j  -  ,  • 

the  Ti.cvpsai’.  num.c  r  i  c  a  1  1  ■■  generated  bodv-1  itteC  la.  rdinati  ■  :  d.,  .  i 

Tilt  iiet'.id  of  attra.  tin,:,  gnja  lint-  te  the  la!  ‘iirfa  ,  a-.'  '.-•i,' 
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surface  points  are  also  discussed.  The  governing  equations  for  incon- 
pressible  two-dimensional  turbulent  flow  in  primitive  variables  are 
presented  in  Section  III.  Boundary  and  initial  conditions  are  discussed. 

An  algebraic  eddy  viscosity  turbulence  model  formulated  for  adverse 
pressure  gradient  separated  flows  is  described.  An  algebraic  type  of 
turbulence  model  has  been  chosen  because  of  the  previous  computational 
successes  for  other  turbulent  flow  problems  and  the  definite  computational 
efficiency  obtained  by  this  approach.  The  method  for  evaluating  forcf 
coefficients  e.xerted  on  the  body  by  the  flow  field  is  then  discussLj.  I'.c 
numerical  techniques  used  in  obtaining  a  solution  are  presented  in  ihctii,:. 

IV.  First,  the  finite  difference  method  for  numerical  I;.'  generatin'  ti.e 
grid  is  given.  The  h'av ier-St okes  equations  are  then  written  in  finite 
difference  form.  The  numerical  rr.etb.od  used  to  implement  the  boundar;.  an.d 
initial  conditions  is  described.  The  computational  procedure  for  introducing 
turbulence  is  next  discussed.  Tlie  finite  difference  and  nutarical  inte¬ 
gration  techniqvies  for  calculating  the-  force  coefficients  are  then  preseir. «.  d . 
The  form.al  di.-.cussion  of  the  numerical  solutions  and  comparisons  v  i  i e:::.:';- 
tiental  data  are  given  in  Section  V.  The  conclusions  derived  fro:-,  tl.i-  resiatn'' 
work  and  recommendations  tor  future  work  art  set  fort!,  in  Sectitm  Vi. 


Partial  derivatives  are  transformed  using  the  chain  rule  and  inver,- 


transformation  relations.  For  a  sufficiently  dif f ei en t iaf. le  function  f 
of  X  and  y,  the  derivatives  transform  as  follows: 

f  =  (v  f  .  -  V  .  f  1/J  * 

X  '  r  i  *  ;  ’ 

f  r  +  ^  '*/•’  ( 

Higher  order  derivatives  are  similarly  derived. 

The  system  of  elliptic  coordinate  generating  equations  is  chosen  to 
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and  on  curve  C  . 
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where  ’(.(x,  y)  are  prescribed  functions  and  and  art 

prescribed  constants.  The  other  two  sides  and  C.,  in  the  rectangular 

region  are  transformed  f rot.  the  branch  cut  -  C^.  The  function.s  x(.',  '1 

and  y(  )  and  all  derivatives  are  continuous  across  this  cut.  These 
boundary  conditions  insure  that  the  body  and  the  far  field  boundaries  ari 
defined  by  constant  ’  lines.  The  generalized  functions  P  and  are  ntili.”, 
in  attracting,  the  coordinate  lines  in  various  regions. 

The  rectangular  grid  system  in  the  transformed  plane  is  used  1 o !'  the 


compu  t  a  t  i  on  s  .  The  mesh  gt'nerating  system,  V.quations  ha  and  ,  is  tratis- 
forried  to  the  co!ti['u t a t  i ona  1  plane  and  becomes 
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TFa-  '  line  spaeini-'  is  adjusted  through,  the  use  of  the  0  funetiini  proj'est. 
!'y  ’lionp-en  ('-1  and  modified  by  Hodge  (Tl 
K( 

0  (  '  ,  ■  1  =  -  i  A,  (  ''  exp'-F(  h  ■  -  ■  ,  1  vl 

k=l 

where  are  tiie  arrp  1  i  f  i  ca  t  i  on  faetors,  I'  is  t’le  dar.piny  fattt'r,  h  is  t', 
number  of  ter-'s,  and  =  (k  -  11  '  .  Tin.  ampl  i  f  i  ca  t  i  t'n.  ar.d  da"  pin.:  tr.cl 
reqiiiring  tliat  the  t  ange::  t  i  a  1  vt  loci  tv  n.  i-e  a  liiiear 
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are  determined  bs’ 


The  second  part  of  the  boundary  layer  grid  line  attraction  techni^ue 


involves  calculating  values  of  and  D  which  when  used  in  0  ^'ive  tl^c 
desired  n  line  spacing.  Near  the  body  surface,  the  transf orniat icn 

equation,  Equation  6b,  can  be  written 

(h  )  =  Q(  ■,  r)  on 

where  variations  with  respect  to  x  are  neglected.  If  the  expression  for 
the  forcing  function  Q  in  Equation  15  is  substituted  into  Equation  20, 
then  Equation  20  can  be  integrated  with  respect  to  to  obtain 
2 

k=l  2(Aj,/I')  expI-D  H('',  +  H(sgn(''-j^  -  0.) 

[1  -  exp(-p'’'  -  )]:  (21) 

where  the  difference  function  H  and  sgn  function  are  defined  as  follows 
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Now  evaluate  the  expression  for  (d^./dy)"  ,  Equation  21,  at  '^  =  0  and 

r,  =  “  ’.■A’'  which  gives  the  approximate  expressions  written  in  Equation 

22  and  23,  respectively. 

'><d7>  0,.o  -  D  \  <") 


dy  .  "  “3- 

may. 


•  1  '' 
k=l 


(23) 


The  derivative  in  Equation  22  can  be  evaluated  using  the  v  values  obtained 


previously.  The  derivative  in  Equation  23  can  be  specified  by  choosing  an 


m 


’“i  based  on  tlie  placement  of  the  outer  boundary.  Then  Equations  22 
max  '  ' 

and  23  can  be  combined  to  give  the  following  equation  for  the  damping 
factor  D: 


D  =  In  2(^) 
dv 


('— ) 
=  0  Mv 


■  .  7.'  -  V-  - 

-  ’  max  K 


(24) 


With  the  damping  factor  D  known,  the  analytic  expression  (Equation  21) 
can  be  evaluated  on  each  interval  of  successive  y  values  obtained  from 
the  Blasius  solution  along  with  the  last  2.!  interval  near  the  outer  boun¬ 
dary.  This  procedure  gives  a  system  of  K  equations  for  the  Y.  unknowns 
.4^  which  are  then  solved.  Therefore,  values  for  the  amplification  factors 
and  damping  factor  are  obtained  at  each  f  line  which  are  used  to  compute 
a  body  fitted  grid  system. 

C .  Body  Surface  Point  Distribution 

The  body  surface  can  be  defined  by  either  a  given  set  of  points  or 
by  an  analytical  expression.  The  body  surface  grid  points  can  then  be 
determined  by  either  interpolation  or  evaluation  of  the  analytical  expres¬ 
sion. 

The  WACA  0012  airfoil  section  is  used  in  this  investigation.  This 
airfoil  has  been  used  in  several  investigations  (65)  and  is  an  AGART  (66) 
designated  test  airfoil  section.  The  KACA  four  digit  series  airfoils 
have  an  analytic  description  of  the  body  surface  given  by  (67) 

y^(K)  =  5f,'0.2969x  ‘  -  0.126x  -  0.3516.x'  +  0.2843x^  -  O.lOlSx'"';  (25) 

where  x  is  the  distance  along  the  chord  (nondimensionalized  by  chord),  y^ 
is  the  upper  surface  ordinate,  and  t  is  the  thickness  given  as  a  fraction 
of  the  chord.  For  the  WACA  0012  airfoil,  t  =  0.12.  The  thickness  function 
given  by  Equation  25  does  not  close  the  body  at  the  trailing  edge  (x=]). 
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A  circular  arc  is  used  to  close  the  body.  The  arc  is  constructed  to  have 
a  center  on  the  chord  line,  to  have  points  of  tangency  with  the  upper  and 
lower  airfoil  surfaces  near  x  =  0.99,  and  to  intersect  the  chord  line  at 
X  =  1.0. 

The  distribution  of  grid  points  on  the  airfoil  surface  is  based  on 
both  resolving  streamwise  flow  field  gradients  and  defining  the  surface 
curvature.  These  criteria  require  the  clustering  of  points  in  the  nose 
and  trailing  edge  regions  of  the  airfoil.  The  following  analytic  trans¬ 
formation  is  used  to  obtain  the  desired  clustered  distribution  of  x 
abscissa  values 

A3 

tanh(— -  tanh(- 

^  =  ;  - A, - £1^  (g.) 

-  Al,  _  Al, 

tanh( — — )  -  tanh(-  — ) 

where  z  represents  the  input  set  of  equally  incremented  values  (0  ^  z  ^  1) , 
X  represents  the  clustered  set  of  abscissae  (0  <  x  <  1),  and  Al,  A2,  and  A3 
are  arbitrary  constants.  The  constant  Al  determines  the  center  of  the 
hyperbolic  function  and  is  used  to  cluster  points  toward  one  end.  The  con¬ 
stant  A2  varies  the  slope  of  the  function  at  the  center  and  thus  determines 
the  distribution  near  mid  chord.  Tlie  exponent  A3  can  also  affect  the 
clustering  of  points  at  the  ends.  If  A3  >  1,  more  points  occur  near  x  =  0; 
while  if  A3  <  1,  points  are  clustered  more  at  x  =  1. 
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SECTION  III 


GOVERNING  EQUATIONS 


The  tirnt.'  dependent  Infompress ible  viscous  Navier-Stokes  equations 
for  two-dimensional  flows  are  presented.  The  primitive  variables  of  fluid 
velocity  and  pressure  are  used.  The  Reynolds  averaged  Navier-Stokes  equa¬ 
tions  are  obtained.  Boundary  and  initial  conditions  are  given.  A  modified 
turbulence  modoi  based  on  Prandtl’s  mixing  length  theory  is  described  which 
expresses  the  Reynolds  stresses  in  terms  of  mean  flow  quantities  for 
separated  adverse  pressure  gradient  flows.  An  expression  for  the  force 
exerted  by  the  flow  field  on  the  body  is  derived.  The  equations  and 
boundary  conditions  are  transformed  by  the  curvilinear  transformation 
discussed  in  Section  II  and  written  in  the  computational  plane. 


A .  Basic  Equat ions  of  Motion 

The  time  dependent  incompressible  Navier-Stokes  equations,  w’hich 
describe  conservation  of  linear  momentum,  in  orthoncrmal  indicial  notation 
are  given  by 


^ -  +  ^.(V.V.).'  =  — 5.p  +  3.  T.. 

J  J  1  1*^  j  ji 


(27) 


where  is  the  fluid  density,  p  is  the  fluid  static  pressure,  are  the 
components  of  fluid  velocity,  and  is  the  viscous  stress  tensor.  The 

subscript  i  =  1  corresponds  to  the  x  dirction  and  i  =  2  to  the  y  direction. 
Conservation  of  mass  for  an  incompressible  flow  is  given  by 


?.  V.  =  0  (281 

2  2 


be 


Turbulent  flow  is  introduced  by  assuming  that  the  flow  quantities 
described  in  terms  of  mean  and  fluctuating  quantities  in  the  form 


V  .  =  V  .  +  V  . 

Ill 


p  =  p  +  p' 


X  .  .  =  T  .  .  +  •  '  .  . 


can 
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dro[>  the  viscous  divergence  tern,  and  kei’p  the  convective  divercip.ci  terr. 


as  a  correction  term.  Thus,  the  appropriate  set  of  equations  becomes 


^  +  r*.  (V.V'  )  =  -  3.  p  +  O.3.V.) 

j  J  i  i  Re  "  j  j  i^ 


D  =  3  .V. 
3  3 


(38) 

(39) 


The  calculation  of  pressure  in  the  primitive  variables  formulation 
must  be  carefully  handled.  The  pressure  does  not  appear  in  the  continuity 
equation  nor  does  a  time  derivative  of  pressure  occur.  Two  approaches 
have  been  extensively  used  to  overcome  this  difficulty.  Both  techniques 
compute  the  pressure  by  utilizing  the  velocity  divergence  D  as  a  correction 
factor.  The  calculated  pressure  is  made  to  satisfy  continuity.  The  first 
method,  introduced  by  C.horin  (6S),  computes  pressure  using  an  iteration 


procedure  with  D  as  the  correcting  term.  The  expression  is 
(s+1)  (s) 


=  P 


;  0 


(40) 


where  s  denotes  the  iteration  number  and  !  is  an  acceleration  parameter 
which  can  be  either  a  constant  or  can  be  related  to  a  successive-over¬ 
relaxation  (SOR)  solution  of  a  Poisson  equation  for  pressure  a.-i  sliown  by 
Hodge  (53).  The  second  method  uses  a  Poisson  equation  for  pressure  derived 
from  the  divergence  of  the  momentum  equation  expressed  in  Kqu;uion  38, 

The  form  of  ttie  equation  is  shown  in  Equation  41. 


D  =  V  •  R  - 
t  — 


(41) 


where  R  represents  the  convective  and  viscous  terms.  if  the  K  expression 

is  algebraically  simplified  with  D  used  wliercver  possible,  thi  following 

simplified  incompressible  Poisson  pressure  equation  is  obtained 

a  a 

r>  =  -  ,  (  nl*)  +  (vP)  +  u*"  +  2v  u  +  v^  +  up  +  vP 
t  X  V  x  X  v  V  X  \ 


; —  (n  +  D  )  -  (p  +  p  ) 

<e  XX  yv  ^xx  ' y v 


+ 


This  equation  contains  both  spatial  and  time  derivatives  of  the  veKuit; 
divergence  as  well  as  the  fluid  velocities  and  pressure.  Again,  continuit 
is  used  as  a  correcting  factor  for  calculating  pi.ssure.  Further  sir.plifi 
cation  can  occur  if  small  variations  in  the  velocity  divergence  are  assu::;. 
In  this  case,  the  spatial  derivatives  of  D  are  set  to  zero  in  Fcuatior:  if 
and  the  result  becomes 

■)  T 

D  +  (u^  +  2v  u  +  v^)  =  -  (p  +  p  )  (■•'i 

t  X  X  y  y  XX  '^yy 

Hodge  (53)  used  both  Equations  42  and  43  and  found  no  apparent  difference 
in  the  flov  field  solutions  for  his  case  of  laminar  flow.  Also,  Chorin's 
method  of  Equation  40  is  related  to  the  pressure  Equation  43.  Hodge  (53i 
showed  that  the  Chorin  technique  is  approximately  related  to  a  solution 
of  the  Poisson  pressure  equation  using  SOR  iteration  on  a  coarse  grid  of 
twice  the  spacing.  The  Poisson  pressure  equation  technique  is  cliosL-n 
for  the  interior  flow  field  because  of  the  Increased  coupling  whici’.  occur.'- 
in  a  finite  difference  representation.  The  body  surface  pressure,  how¬ 
ever,  is  calculated  using  the  iteration  technique  of  liquation  40.  Hi-rL 
mass  conservation  is  imposed  directly  and  the  difficulty  of  evaluating 
the  Laplacian  of  pressure  at  the  bod;  surface  in  Equation  43  is  avoid^■d. . 

The  following  set  of  equations  taken  from  Equation,'^  38,  3'^  ar.d.  4a 
becomes  the  system  used  for  the  solution  of  inconipress  ibl  e  iurbulei;t 
viscous  two-d imens iona 1  flow  over  airfoils: 


u 

t 

4. 

uu  -f  vu  +  uP  = 

X  y 

1 

u  ) 
vv 

Re 

t 

vu 

XX 

{  •  4  •  .  ' 

V 

t 

+ 

uv  +  vv  +  vD  = 

X  y 

-p,  + 

1 

Ke 

t 

(v 

XX 

+ 

V  ) 

V  V 

D 

t 

+ 

-f  2v  u  +  = 

X  XV  V 

.  +  1 

1  ) 
vv 

n  =  u  +  V  ( .t  7  I 

X  \’ 
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nu::if  r  i  c  a  1  solatiun  ol  t  Vu-  L’OVi-rninj.-  t-qu.i  l  i , 'ii 


forr.ing  the  conpu  t  a  t  i  ans  in  the  eo::.put  at  ional  iilaiii  nl  tainih  vit’.  t 
vilinear  transformation  discussej  in  Sent  ion  '1.  !:.i^  •eth.d  ri's.ir 

the  equations  are  transformed  to  the  comput  a  t  i  i.'na pjani  .  fhi  »-L'/er 
equations,  ?!qnations  44  tfirono.h  47,  art  t  rans !  or."-ed  u-ir.o  the  r<  ’a;i 
in  Section  II  and  Apjiendix  A.  The  resuitinj;  set  <  f  c  r  it  -  p  J  : tr 
formed  equations  hecon.es  : 


U  ,  V  c  / 

U.J.  +  -j-  (y,  u  f  -  y  .  u  )  +  q  ( X  .  u_ 


1  /  '  ’  ( 
J  p  :  -  y :  i>. »  + 


+  u,  (---)  +  u  .1  - , ; 

d"  '  ,r 


+  —  (>q  V  f  -  y  -  v_.  )  +  '-  ( :•;  .  v_  -  v  .)  + 


-  -  (X  ;  p,  -  X.  r  d  + 


1  ( 


+  V,  i  +  V  A  -  ,) 

.r  ,r 


1)  I  +  •  ,  ( u  .  -  y  .  u  (  -.q  V 


4  (  X  .  V  -  X  V  . '  ■ 


P.  <  -t'  +  P  -1 

r  .1^ 


=  (■ 


+  (:•; 


i-. .  Boui'ida  rv  and  I II  i  t  i  a  1_  C.ond  i  t  K'li^ 

Tin’  complete  definition  of  the  flow  problem  requires  the  t-pe^  i :  i^atii  : 
of  sufficient  boundary  and  initial  conditions.  Tlie  ‘.'av  ie  r-S  t  ekes  equation, 
are  coupled,  nonlinear,  first  order  in  time,  second  order  in  spa,e,  illipt  i 
partial  differential  eijuations.  Ttie  elliptic  nature  I'f  the  ecu.ition-  for 
v'eiocit  .'  requir>.‘s  biHindary  conditions  at  ea..h  boimd.ir'.  .  'he  first  orahr 
thricatives  of  pressure  require  a  boundar;.-  condition  descrll  in/  the  fria- 
slrear.  pressure.  Tiu'  time  derivative  terms  require  inili.jl  cone,  i  t  i  on  on 
'.'elo,  it'.'  and  pressure  throughout  t!ie  flow  tield. 


'unw..f'.  c; 


,'r.ditions  on  tin  airfoil  surface  ari-  i-on-  i>h  •  ir  a  . 


eon, :  inuu:-.  n,;  slip  boundar'.'  condition  fi'r  a  viscous  fluid;  at  a  station...;" 
solid  sur:,ioe  is  in;po.sed  oi;  t  .he  tangential  Velocity  at  t  !u  surf.ue.  1;' 

a.i.l  ;  t  i  ,  the  nor;'.al  component  of  velocity  must  he  zvre  ''or  a  st.it  r 
surface  with  no  t  r.iiisp  i  rat  ion  .  These  conditions  are  ’'.et  by  tin  ,s;.ei' i  f  icat  i 
c>:  t  i.t'  biriti.Ut  biHindar'.'  conditi.'U'^ 


v,  e  (. 


w'.L  to  tin  su'  ooi'ipt  b  den.  te^  iiou'.  sul'f.ne.  i.'.e  prissUfi  on  li.i  .a  ;  )  1  .' i  1 
surt.ice  i .-  ^'or.puted,  a:d  n.  bo.'b.'  lunnuiary  conditio:,  is  ft  .uin,!  'n',  t!,..  tie,- 

ii(.  r  ; '.'.1 1  i '.'o  pre~-uri  ter:  s. 

'h.  tar  fiild  boundary  con.!  i  t  i  o::-,  for  VLl.'cit;.  .ind  jrLo.nrt  spi  .  ;  : 
till  free’  trea-,  f  !  o'.-;  ;ii  Id.  let  r  hi  .i  po.'iti.'U  v'e.t.T  in  t'u  fl.'W  ;' i .  Id 
wit'  til.  origin  a'  tiie  .lirfoil  n.id-ihor'.  Ih.-n,  ,i-  t  h.  dist.i'ui.  r  !  r.’--,  ti 


a  1  I'  1  1  1  n 


1  a  r  .u.  t  hi'  fluid  '.'.lo.  ity  .ind  ;  risiiiti-  .ij-pri'ach  t  h.  ir 


[uuiiii-cno  ion.i!  fre-  --trea:'  '.'aliie 


.  i  -t  s  in  . 


w'.i  ri  ir  ti.i  cl  dlta.i-  i  hurti  wi;^  riapi-t  ti 

;  rt'i  t.  ri  a. . . 

:!a  ilr-t  aiji  r  !  di  r  i  va  t  ivi-  ttrKs  recuin-  initial  i  oiui  i  t  i .  T:  -  r^i 
Iht  vr!  -  it',  ami  ;>r>^^\.ia  in  tat.  field.  Fra:,  a  tlaariliial  viLVp.  i::t, 

--  i  'a  a  t  ■  I  a '  ;  t  t  i  a  ~  t  ,  ad;  1  u '  i  un  is  des  i  Tt  si  ,  a:  a.  ma  t  i  nn. '  ;  i  :  r  ;  a  1  s  ;  i  . 

aatia:  a:  Val.ait'-  a::d  pra-.suri-  (suF;eal  ta  iFa  p  r  a  sc  r  i  la-d  •'iM:;dar;  ^anditaa.  ’ 
is  i a  r".  i  s  1 :  1 1  . 

I  Fa  bL'undar'  tondit  ion.-  and  initial  laaid  i  t  i  rais  ari  .ils^  t  ra:;  ;  ar::a  a’  ta 
t  Ff  ca:-:;'n  t  a  t  ; '.'na  !  pla-u.  iFa  valuas  lor  aa^i,  aa:-;.;;;;'’.  rt-'rair  \:Tialttrad  an: 
an  tFa-  :a  1  rapi-atad  lati.'.  'i:.a  loiatio;;  i:;  t  Ft  :  ra:n  :  a  r:,.f  J  ;a,ai  dt;a::dt  n: 
t:if  -pf.  i:iv  :  rans  :  a  ta.a  t  i  o;:  t  ha^ta.  to  pt-nt-ratt  t!,«.-  ooraa ;  t  a  t  i  ca:..  1  p  r  i  d  .  dFi; 
r.attt  r  i  d  i  sau.^st  in  Stt'tian  1'.'  wF<  re  t'.e  nu::.tTit,il  I'rot’t  dure.'  art  dfstri’td. 
'.he  stais  i  t  iv  i  t an.ilysis  t'f  the  euter  ht'undar;,'  laa.atiar,  is  a  dd  ressta!  in 
;\ppen..l  i  . 

r.  lurFuK aua;  ;nuK_l 

!  iirFul  t  nae  rsuit.' 1 1  i  i;  ji  is  t'i  considerahle  iriti-rest  in  theoretical,  ey.]u  r  i - 
rantal,  and  oor.putat  it'na  1  research,  wors.  ‘d.iir.-  ,it  rt.>dy::a::.ic  prtblta:,s  v'f 
pra.tia.il  interist  ciM'.t.’.in  rc'pior.s  of  tvirbultait  flow.  Several  r.tsltl'  '..ivi. 
beta:  deve  1  tiped  whicli  varv  in  cor.pl  ex  i  t fro-  I’randtl 's  lta:.t''  t’.i.r 

ft'  t  (W  F.n  i  (jues  employing  au.xillary  d  if  f  ert  a:t  ia]  ee;<-:ation  for  It:;.:' 

visca.sitv  coefficients,  and  turbulent  kinetic  entreits.  I'Fe  atuur.ik  ;  of 
each  riodel  varies  with  both  the  geor.etry  an!  t ’n,  f  K'w  ;i,  ;d  ’  ;.i: .  l  i  I  ■.  a,.,'  l,  t 
Ct'r.p.iring  results.  in  computational  work,  tht  .ibilit'  te  tr.i:,  i:.:.  :  F,  :i:i- 

bulencc-  mitdel  into  numerical  procedures  and  tlie  taa-puter  re- -unt  n  ai;:,.! 
are  aiidition.al  factors  which  need  to  be  con.- i  dett  il . 

In  this  inve.st  igat  ion ,  turbulence  is  modtlltt!  with  a:i  al.u'r..,.  t 
visct'sity  approacb  based  on  Prandtl's  miy.in,-  le:;gth  thn'r-  .  ui  '  ;• 

concept  of  tliis  mt  thod  is  tiiat  tht  turbulent  Kt-.noKl  at  ru  ,  a  ,  -  . 


L-ar.  hi'  rt-lati'd  ti>  tln'  strain  rate  of  the  mean  flov  in  a  va-.-  anala^’otis  to 
tliat  for  lai’iinar  viscous  stresses.  The  f unda:;;en ta  1  expression  was  previous);, 
introduced  in  Kcjuation  jd  to  obtain  the  simplified  turbulent  incor.press  i )' 1 1 
Na  V  ii' r-S  t  okes  etjuation,  Equation  13.  Tlie  algebraic  eddv  viscositv  lui.- 
been  used  exti'iisivelv  in  obtaining  numerical  solutions  for  a  varief  of 
aerodynar..  i  c  pro)  lerr.s  (  d  A  ,  d  3  ,  *'•9 1  . 

Dll'  edd;.-  viscositv  model  for  turbulent  flow  wit'i  r'ero  pressuri  'cradieut 
near  a  Si’lid  surface,  which  incorporates  the  modifications  of  ('eb.c  i  and 
S::'.it)i  it-ii.  Van  I'riest  (di),  and  Pelwert  (37'),  is  given  below  and  then  discus 


(k,  vl'  >  (u"  +  v"l 

11.  X 


=  1  -  exp 


1  +  3.3  (y/D 


(1  -  ^-1  dv 
u 


The  tancential  surtaci'  diri'ction  is  given  by  x  and  t  hi'  normal  s’arlaci 

direction  le.-  v  .  rile  inner  and  outer  I'dd'.'  v  i  scos  i  t  ii's  .  ■  and  .  modi  ) 

i  c 

till  law  of  tile  Wall  and  l.aw  of  the  U’ake  (7('’)  recions.  Tin-  inner  mixii;.- 
leru-th  is  assumed  proportional  to  the  distance  v  from  tlie  surface. 
the  Von  Karman  constant.  The  Van  Priest  (dD  dampinc  factor  r  oka  1  ■ 
till  laminar  suMaver  region  when-  is  the  surfai  i  slu-ar  stress  and; 
the  decay  ci'iistant.  The  displacement  tliickiii-ss  '  iietiiud  bs  Enuation  ' 


is  the  outer 

la; 

■er  len, 

■th  scale.  I  he 

i  ni  t-r:'.  ill* 

liu-y  ‘.iv'i-'r  ■, 

ih  ; 

in  Equation 

nodi'  1  s 

tlu'  decreasi-  of 

t  :u'  tiiris 

s  K‘:k' t  •  i  :  iT.  -  i 

t  c  f 

fully  turbul 

t‘I1  t 

region 

t  o  t  lie  i  nv  i  sc  id 

rt’.  ie’n  .j; 

^■ei^  {}■,  !> 

The  intermit 

,y  is  a 

f  line  t  i  on  (1 :  ;  hi- 

nurr.<ii  f! 

i  .j:d  (■  t  Ti-r  : 

: ,  1  Vs 

nond  ini-ns  ional  iat'd  by  l)ie  local  boundary  layer  distance  Ii.l- 

edi'e  boundarv  laver  volocif.  u  is  the  outer  laver  characteristic 

-  t. 

Ihc  standard  (  cl'ec  i-Sitith  values  for  the  constants  are  =  .-11, 

k,  =  .(M*  a:;d  A  =  2<:.  Continuity  of  the  eddy  viscosity  d  i  s  t  r  ibv:  t  ;  o'c. 

achiived  b;:  sv  it  chine  iron,  the  inner  to  the  outer  r.odel  vb.(_ne\Lr  ,  fir 
esoevd- 

K._'Centlv,  Jo'  e  and  Hankey  (71)  compared  this  eddv  viscositv  rtc'de; 
b-.  Ihiuat  ions  55  throuch  S'J  with  experimental  data  for  several  constant 
adverse  prL-ssiire  gradie>!t  attacined  flouts.  Tisev  varied  t!ie  k. ,  ,  k,,  ai.d. 
parar,eti.rs  in  a  numeric  il  turbulent  boundary  layer  prograr,  until  s, 

tion.s  matched  t'^e  expe  r  incnital  boundary  layer  velocity  profiles  and  ski 
friction  cOe  f  ic  ien  t  s  .  These  pa’-ameters  had  substantially  different  va 
than  tlu  previously  cited  values  obtained  for  zero  pressure  gradient  fl 
They  found  ctiat  when  the  flow  encounters  an  adverse  pressure  gradient  t 
eddy  vi.sc■Osit^'  para.-neters  rapidly  adjust  to  the  values  of  =  .b5  and 

4. 

A  =  3d.  The  outer  eddy  viscosity  constan.t  appears  to  linearly  increas 
wit;,  distance  measured  from  the  start  of  tbie  adverse  pressure  gradiei'.t. 
Tluis  .Iob'i.‘  and  hankey  (71)  found  that  k,,  can  vary  from  a  value  of  .Clf' 
valut-s  of  .1)3  or  greater  given  that  the  adverse  pressure  gradient  e:-;ist 
for  a  sufficient  length  downstream. 

Ihi;  use  of  an  edd^■  viscosity  model  similar  in  form  to  the  on<-  givt 
in  hjuations  55  throv:gb.  39  for  an  airfoil  at  angle  of  attack  recuirt,- 
examinat  ion  of  the  topical  fl.'w  field  characteristics.  Tb.t  upper  sr.rfa 
of  the  airfoil  has  a  pressiire  distribution  wb.ich  includes  a  regioi’.  of  1 
pressvire,  the  sv.ction  pea.k,  near  the  airfoil  lc>ading'  edge.  The  .sin  •  i  s. 
is  followed  bv  an  ad'c.  r:  i  pressure  gradient  ri.clL'n.  Tin.-  pri  s.-^o.in.  ,  i  .A  i 

lend-  to  dt-  rea,-e  a  thie  ‘"low  recovers  toward  tbie  trailing  (v’.gi  . 
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turbulence  nociel  that  is  used  for  an  upper  surfaee  turbuli  nl  b  jin.c.c:' 
should  include  the  effects  of  this  adverse  pressure  pradier.t.  i!’, 

eddy  viscosity  parameters,  therefore,  take  C'n  the  Johe-  ar.c  iinnl.e  i7] 

of  =  .bS  and  a"*”  =  52  in  this  region.  l!ie  outer  odd"  visce.'if  ]  .ir 
k.-i  ntav  alsv)  change  in  th.e  adverse  pressure  gradient  region.  i;  ■'...c'Vc  r  , 
rate  of  chance  is  affected  bv  the  pressure  gradient  ar.d  dovnstrear.  di 
The  pressure  gradient  does  decrease  rapidly  as  the  flov  proceeds  dovn 
.Any  tendeiicv  for  k>  to  increase  is  probably  offset  hy  ti-.e  decreasir.c 
pressure  gradient.  Ihc.-,  tht‘  value  of  k,  is  kept  at  thi'  ecnili:  r  i  a:t 
of  .  l !  1  '■  b  . 

The  ei'dy  visc.sity  m.sh  1  wit!,  the  discussed  adverse  pressure  era 
modif  icat  ion.s  i-  now  written  in  terr..-  of  t'.e  nond  ime!’’..-;  ional  varial'les 
duced  in  r.i'p.ni  t  i.jn.s  f-<  and  37.  .All  length  epaantities  arc  nend  ira.  n.s  ior. 
with  respect  to  the  airfoil  ciiord.  'iiie  mea"!  and  nond  imens  iona  1  quant 
are  hencefortli  implied.  Tlie  eddy  viscosity  model  I'ecor.es 
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c.  anc;  are  considered  adver.se  pressure  gradient  paratete 
de'v  viscositv  node]  niist  a].s[>  r ejire.seit t  t!u  phvsicnl  ciiarac 


ii.  l.irce-  se'csi  r,i  •  ed  re.:  ions 


Tra  ilia 


which  e.ccur  during  .airfoil  stall. 


stall  occurs  when  the  turbulent  boundary  layer  separation  r;oves  tovard  the 
leading  edge  of  the  airfoil  as  the  angle  of  attack  increases.  Ihe  pressure 
distribution  in  a  trailing  edge  separated  region  has  snail  gradients  (72). 

This  observation  motivates  another  changi-  to  the  eddy  viscosity  model.  The 
inner  eddy  viscosity  parameter  is  relaxed  hack  toward  the  equilibrium  value 
of  .41.  The  separation  point,  which  precedes  the  return  to  zero  jiressurc  grad¬ 
ient,  is  selected  as  the  point  to  start  the  exponential  decay.  Ti.e  following 
relaxation  form  with  distances  nondinensionalized  by  the  chord  is  used 

s  -  s  , 

f  =  .'i  1  +  B  exp(-  - )  s  •  s 

r 

I 

f  =  1  s  ■'  s  . 

d 

where  f  is  the  relaxation  factor  which  multiplies  the  initial  value  of  kj, 

k“^,  s  is  the  distance  downstream  from,  the  trailing  edge  separation  point, 

s,  is  the  delav  distance,  s  is  the  relaxation  distance,  and  A  and  B  are 
d  •  r 

constants.  The  constants  A  and  3  are  determined  by  applying  the  conditions 

2 

f  =  1  at  s  =  s^  and  f  =  s  approaches  infinity  where  is  t}';e 

asymptotic  value  of  dow’nstream. .  Thus,  A  and  B  becomc- 

A  =  B  =  A~^  -  1  ('^6) 

Measurements  of  Bachalo  (73)  and  Baker  (74)  in  the  trailing  edge  separa¬ 
tion  region  indicate  the  presence  of  a  lower  but  finite  turbulent  inters  itv 
when  compared  with  the  separated  turbulent  shear  layer.  Tlie  velocity  deriv¬ 
ative  terns  in  the  inner  eddy  viscosity  and  damping  factor  expressions,  Equa¬ 
tions  60  and  62,  can  approach  zero  in  a  region  near  a  velocity  inflection 
point.  In  this  case,  a  zero  eddy  viscosity  is  perm.itted  in  a  turbulent 
separated  region  which  is  contrary  to  the  apparent  flow  structure.  T!iis  same 
zero  eddy  viscosity  behavior  can  also  occur  in  the  case  of  .cheddine  .separation 
bubbles  from  the  airfoil  surface.  In  order  to  prevent  tlie  appearance  of 
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laminar  regions  within  the  turbulent,  adverse  pressure  gradient,  separated 
flow  field  during  airfoil  stall,  a  limiting  technique  is  employed.  The  inner 
eddy  viscosity  is  prevented  from  decreasing  in  value  as  the  normal  distance 
from,  the  surface  Increases.  This  restriction  simulates  the  uniform  turbulent 
intensities  measured  (73,74)  in  the  trailing  edge  separated  region.  Tlie  innt;  !' 
eddy  viscosity  is  further  limited  by  imposing  a  condition  of  no  decrease 
in  the  do\,mstrean;  direction.  This  limit  provides  for  a  finite  turbulent 
intensity  in  separated  regions  which  may  develop  in  the  turbulent  boundary 
layer.  The  limit  is  initiated  by  obtaining  the  distribution  of  eddy  vis¬ 
cosity  in  the  attached  boundary  layer  near  the  leading  edge  of  the  airfoil. 

This  distribution  is  compared  with  the  distributions  calculated  by  the 
model  at  downstream  locations.  The  outer  eddy  viscosity,  which  has  no  velocity 
derivative,  is  not  limited  in  any  way.  The  "limiting"  technique  is  passive  in 
the  sense  that  the  locally  computed  value  for  the  inner  eddy  viscosity  is  used 
whenever  possible.  Numerical  experiments  without  the  "limiting"  technique 
displayed  unphysical  results  because  the  conventionally  modelled  eddy  viscosity 
is  divergent.  The  inception  of  separated  flow  decreased  the  eddy  viscosity 
thereby  causing  futher  separation.  The  "limiting"  concept  is  necessary  to 
prevent  this  unphysical  result. 

The  transition  from  laminar  flow  to  fully  turbulent  flow  is  modelled 
using  the  transition  model  of  Dhawan  and  Narasimha  (75).  The  expression  for 
the  transition  factor  T  is  given  below 
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model  are  conveniently  expressed  in  terms  of  quantities  in  the  transforriC-d 
plane  described  in  Section  II.  The  velocity  derivatives  in  tlie  inner  eddy 
viscosity  model  are  formulated  by  determinii.g  the  tangential  and  normal 
components  of  velocity,  u  and  v  respectively,  to  the  body  surface  ^ 
at  a  position  ' 

u(  r,  =  c.,u(  r,  r)  -  c^v(>,  •  )  (bb; 

v(  r)  =  Cj^u(  c,  y)  +  c.jV(  5,  r,)  (69) 

where  and  c.,  are  components  of  the  unit  outward  normal  to  the  body  sur¬ 
face  at  position  (  given  by 


c,  =  -V  ,/y  ,  I  — 

Then,  the  directional  derivative  definition  is  used  to  obtain  the  surface 
normal  derivative  of  u  and  the  surface  tangential  derivative  of  v  expressed  as 


-■  u  •  , 


-  V  -U 


. )  +  c  -  x,.u  -)  /J 


(70) 


=  ■c,(y  V,  -  y  ;V  )  -  c,  (x  o',  -  X  v.)./J  (71) 

which  correspond  to  the  u  and  v  derivatives  of  the  Cartesian  eddv  vis- 

^  y  X 

cosity  model  expressions.  Equations  60  through  64. 

The  turbulence  structure  in  the  far  wake  is  modelled  with  an  eddy 
viscosity  expression  similar  to  the  outer  eddy  viscosity  model  in  the  tur¬ 
bulent  boundary  layer.  The  model  is  based  on  the  assumption  of  a  self¬ 
preserving,  equilibrium  turbulent  wake  with  constant  pressure  in  incomprii's- 
sible  flow.  The  approach  has  been  experimentally  investigated  by  Narasimiu. 
and  Prabhu  (76)  who  obtained  the  following  expression 


=  .  k.,w 
w  i  o  w 


(721 


where  w  (x)  =  max  ' b  -  u(x,v)  ■  is  the  maximum  velocity  defect  in  the'  wake, 
o 

V 

is  a  half  wake  width  defined  where  the  velocity  defect  becomes  one  lialf 
w 

the  maximum  value,  and  k,^  is  a  constant  equal  to  .Oeb.  ]f  the  exper  iiiu'T',  t  .i  1 1 ; 
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integrated  value  of  the  wake  defect  function  is  used,  then  Ecjuation  72  can 
be  written  in  terms  of  the  parameters  of  edge  velocity  and  displacement 
thickness  similar  to  boundary  layers  as  follows: 

—  =  k,u  5*  Re  (73) 

3  e  V 

where  u  is  the  nond imensional  velocity  at  the  edge  of  the  waVie ,  :  is  ti.c 
e  w 

nondimensional  displacement  thickness  for  the  half  wake  as  defined  by  Equa¬ 
tion  64,  and  k^  becomes  .0634.  This  expression  has  been  utilized  by  Green, 
et  al  (77)  in  an  integral  method  for  predicting  two-dimensional  incompressible 
and  compressible  turbulent  wakes  of  lifting  airfoils.  For  an  incom.pressib.le , 
constant  pressure  wake  the  momentum,  integral  equation  reduces  to  a  statemient 
that  the  momentum  thickness  is  a  constant.  These  conditions  exist  in  t 
far  wake  of  an  airfoil.  In  addition,  in  the  far  wake  the  displacement  d 
momentum  thicknesses  become  equal  as  showTi  by  the  experimental  r.  .ults  of 
Narasimha  and  Prabhu  (76).  Thus  the  far  wake  eddy  viscosity  is  approximate! v 
constant.  The  transition  in  the  wake  from  the  turbulent  core  regior.  l.o  the 
Inviscid  freestream  region  on  each  side  of  the  wake  is  modelled  with  the 
intermittency  factor  of  Equation  63  where  f  is  now  the  wake  half  width  and 
y  is  the  normal  distance  from  the  wake  center. 

The  turbulent  near  wake  region  mixing  rates  where  the  wall  boundary 
layers  merge  and  eventually  become  the  far  wake  are  not  well  understood  fro:;: 
either  an  analytical  or  experimental  viewpoint.  In  the  flow  regime  undi-r 
consideration,  a  laminar  boundary  layer  from  the  lower  airfoil  surface  r.  i:-us 
with  a  turbulent  boundary  layer  from  the  upper  surface.  An  eddy  viscosity 
model  similar  to  that  used  in  the  outer  turbulent  boundary  layer  and  tlu-  far 
wake  is  used.  Inouye,  ct  al  (78)  applied  the  Cebcc i-Smi tli  (G3)  outer  laver 
model  in  the  near  wake  of  a  flat  plate  and  obtained  good  agreL-an'  for  the 
mean  velocity  profiles.  Tliey  allowed  the  constant  k.,  t('  reach  twice  tb, 


equilibrium  value  of  .0168  over  a  length  of  one  chord  from  the  trailing 
edge.  The  computed  eddy  viscosity  was  still  about  50'  low  when  compared 
with  experiment  at  the  one  chord  distance.  Later,  Burggraf  (79)  compared 
several  turbulence  models  for  various  near  wake  flows  and  concluded  that 
the  Cebeci-Smith  model  adequately  predicted  the  mean  velocity  prefiles. 

The  merging  of  the  upper  and  lower  surface  boundary  layers  is  modelled 
with  a  simplified  interaction  hypothesis  approach  introduced  by  Bradshaw  (80). 
The  two  layers  are  initially  allowed  to  develop  downstream  a  distance 
without  interaction.  At  this  point  the  eddy  viscosity  is  increased  exponen¬ 
tially  in  the  adjacent  laminar  boundary  layer  from  a  zero  value  to  t!,L'  far 

wake  calculated  value  c  at  a  distance  WL„  from  the  trailing  ec'e,  This 

ow  2 

eddy  viscosity  variation  is  formulated  as  follows 

s  -  WL 


c  =  1/2  £  1  +  tanh 

ow 

c  =  0 


(8 


UL.,  -  WL^ 


)■• 


V.T, 


s  1 


(79) 


where  WL  is  the  distance  from  the  trailing  edge  to  a  locatitn  half  way 
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between  the  distances  and  WL.^  and  s  is  the  local  distance  from  the 
trailing  edge.  The  coefficient  8  is  selected  so  that  •.  0  at  s  =  V.'Lj 

and  £  -  at  s  =  because  at  s  =  the  hyperbolic  tangent  function 

becomes  equal  to  0.9993  while  at  s  =  WL^  it  becomes  -  0.9993.  The  eddy 
viscosity  profile  in  the  turbulent  boundary  layer  near  tlie  trailing  ed.-.e 
is  extended  into  the  near  wake. 


D .  Force  Coefficients 

The  force  which  is  exerted  on  the  body  by  the  moving  fluid  arounti  tin' 
body  can  be  determined  from  the  Cauchy  integral  equation  for  consL-rvai  i  ^ni 
of  linear  momentum  applied  at  the  body  surface  or  on  a  closed  ['atli  surround¬ 
ing  the  body.  Tlie  basic  expressions  for  the  force  coefficients  f  orir.u  1  .i  t  ic! 
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,-2x  f-p 


where  the  vorticity  on  the  body  surface  is  given  by 


-  (v  .V  +  X  aj  )  /  J 


The  corresponding  mor.ient  (positive  counterclockwise)  about  a  point  P 

located  within  the  bodv  cross-section  with  coordinates  (x  ,  v  )  is  given 

P  ■  P 


/na 


(x  -  Xp-K-Zx.p  -  -  (y  -  yp)(2y.p  -  — 


-)'d-  (^2) 


Next  consider  again  the  more  general  case  of  a  closed  path  around 
the  body  whose  force  coeftiei^nts  are  given  by  Equations  n.l2a  and  D.ldb. 
Assume  that  the  closed  path  is  an  t  contour.  Tlien ,  using  the  transformed 
expressions  in  Appendix  A  for  the  outward  normal  vector  to  an  ’  contour, 
integrals,  and  derivatives,  the  force  coefficients  for  the  hody  in  term.- 


of  quantities  on  a  constant  '  line  in  the  flow  field  become 


-  V  (x  .V  .)  ]  +  ?u(uy  .  -  vx  .)  d 
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where  n,  •"  r;  <  n  .  The  terms  in  the  line  intecral  represent  in  order 
b  —  p  —  max 

pressure  forces,  viscous  forces,  and  convective  outflows  of  linear  momentu:::. 
The  area  integral  term  is  the  time  rate  of  increase  of  linear  momentum  con¬ 
tained  in  the  control  volume  bounded  bv  the  bodv  t,  contour  and  the  selected 

'  '  b 

ti  contour. 

P 

The  corresponding  lift  and  drag  coefficients  are  then  ca]culatec  usino 
Equations  D.13  and  D.IA. 
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SECTION  IV 

NUMERICAL  NETHOP 

Tlu-  numerical  methods  used  to  obtain  solutions  tor  iiu  o:-.;  n-s  ii  i 
turbulent  flow  over  airfoils  are  presented.  Ihi-  nu:;.eri^a'.  iiniti.  div.ere 
procedure  w-luch  yields  the  grid  transformation  is  first  dis<-U^sLd.  i:.e 
implicit  finite  difference  method  for  obtaining  approximate  solul  i  an,-  tk' 
the  unsteadv,  incompressible,  turbulent  Navier-Stokes  ecu.i  t  i  ons  is  r.i'-l 
described.  liuru-r  i  ca  1  boundary  conditions  and  initial  conditions  are  i  - 
sidered.  The  nur:erical  procedure  for  incorporating  the  turbulent  ei,;.!; 
viscosity  model  is  discussed.  The  computation  of  the  force  coef  f  i  c  i  e  n  t 
on  tlie  bod;.'  is  then  presented. 

A  .  G  r_i_d  J  r  aj^s_fo  rmat  ion 

The  numerical  solution  of  the  governing  equations  for  thi.'  t  rarsfir::  a 

tion  is  carried  out  in  the  transform.et!  or  computational  plane  k'n.  a  scyiare 

mesh  (  •',  ’' )  svster..  The  bodv  contour  C,  and  tin.’  far  field  houndav.  f 

h  ■  1 

become  constant  ’  contours  (Figure  1).  The  f.ar  field  boundary  i  .s  usually 
used  to  approximate  infinity  in  incompressible  flc'v.  proMe-,s.  Con.st.tnt  ' 
contours  in  the  physical  plane  are  required  to  transfor":  to  equi-spaaed  ' 
coordinate  lines  with  unit  step  size  in  tlu  transformed  plane.  The  .-.ptu  i 
contours  in  the  phy.sical  plane  i.s  controlled  by  the  set  of  elliptic  ecu,;; 
The  '  contour  spacing  on  the  body  surface  and  outer  boundar  .-  in  tiit  id.v,-  i 
plane  is  determined  by  the  cho,sen  point  distribution  on  each  ,u\irf,iki  .  Ih 
constant  ’  contours  are  also  required  to  transform  to  equi  i -.spat a  tl  ’  ciu'r- 
dinate  1  ine,s  with  unit  ,ste;'  size  in  the  t  ran,sf  ormed  plane.  Tin  lini  , 
are  denoted  by  the  subscript  i  witli  range  1  ^  i  '  ]M,yy  ^here  IN,','-  i  ,s  the 
of  ■  lines.  Similar;-,,  the  '  lin-s  are  identified  by  a  subscript  i  witli 


1 

.i  1 


JIU 

r.in 


41 


of  1  ■  j  •  JM\>.  whore  JMVf  is  t  tie  nunbi-r  of  '  lin<.;-.  'ho  voliu'  j  =  1 

ar.o  j  =  .IMAX  are  tiio  bodv  and  outer  boundary  surfavo^;,  rospoo  t  i  Vl  !  ,  . 

Tlio  constant  ?  linos  i  =  1  and  i  =  IM."iX  denote  t>io  i'rani ’i;  l  ut  in  I  lio 

piiysical  plane  and  are  therefore  identical. 

T:.e  governing  elliptic  transformation  eeiuations,  I.cuatio:'-  or  ^ 

are  written  with  finite  differences  In  quasi-1  inear  izid  lorv.  ft)r  t  loc¬ 
ation  (i,  j).  Second  order  accurate  central  differences  riven  in  Appi.'a;;v 
are  used  for  evaluating  each  term.  The  highest  order  derivative  lor.tai;.. 
the  unknown  at  location  (i,  j)  while  the  t  rans  f  or::;a  t  ion  coe  f  f  i  c  ie:.  t  s  ,  wh  i  o 
contain  lower  order  derivatives,  are  evaluated  witli  previtn;s  values  of  t ho 
variables  x  and  y.  The  difference  expressions  are  gi'.-sr;  i'tlow  ir’ierc  ti,e 
indices  i  and  j  are  understood  when  omitted. 


X  .  ,  =  ( X  .  ,  ,  +  X  .  , )  -  ( X .  ,  .  ,  ,  -  X  .  .  .  ,  +  .  ,  .  ,  -  x  .  ,  ,  .  ,  1 

ij  1+1  1-1  d  1+1  j+1  1+1  j-1  1-1  j-1  1-1  j+] 
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+  +  v._^)  +  -  y,_,lr  +  /2f  ■  + 


T!ie  i  index  near  the  brancli  cut  is  adjusted  as  folli^ws: 

when  i  =  1  then  i  -  1  =  IM.AX  -  1 

wiien  1  =  IMAl!  -  1  tlien  i  +  1  =  1 

when  i  =  IMiX  then  i  =  1 

Tile  set  of  finite  difference  equations  lias  the  range  of  indices  given  h-, 
1  ■  if'  Ih.'cX  -  1  and  2  j  j_  ,IM-\X  -1  which  results  in  (  1  M'd’- 1 )  (df.A'- -  2 1 


equations  for  the  unknown  values  of  x  and  y  at  the  interior  grici  p’oints.. 


1 


Ihc  syste:-.  of  difference  equations  is  solved  usiny,  succtcsive  over  relaxa¬ 
tion  (SOR)  iteration.  For  a  general  function  f,  an  appro>;in:al  e  solution 
for  f  using  SOR  iteration  has  the  form 


f  =  .f"  +  (1  - 


:  h  / ) 


where  s  represent.s  the  iterate  number,  is  an  accelerat  ion  j->arat;et  er ,  and 
*  represents  the  current  value  of  f  calculated  from  the  difference  equation. 

The  acceleration  parameter  is  always  one  for  the  first  iteration.  The  syster. 
t'f  linearisid  difference  equations  giv«_n  by  liquations  S5  and  86  obtained 
fro:-,  the  elliptic  generatinv  differential  equations  is  consistently  t'rdered  (  ■  1 
i'ius,  loc, -il  optinu:-;  acce  I  erti  t  ion  parameters  ...  can  be  calculated  fro::. 

tj 

theory  (81)  and  are  given  by  llodgt  (nj).  Convergence  of  the  iteration  procedr.ri 
is  determined  by  com;parinr,  the  ab.solute  value  of  the  difference  in  succe.'-'-- i\'e 
iterates  witli  a  prescribed  error  criteria  T. 

SOR  iteration  technique  requires  an  initial  guess  for  tiie  solution. 
Geometric  contours  wit:;  .i  shape  similar  to  the  outer  computational  hound.ir;.- 


arc  used  for  constant  ’  lines  .and  .straight  lines  which  enan;ite  f rot.,  tb.t  bed; 
center  arc  used  for  '  lines.  .-\  similar  appro.ich  has  been  successfully  used 
by  Thames  (4d)  and  Hodge  (63). 

The  forcing  function  0(’,  ’)  defined  by  liquation  13  is  l.•valuatl,  d. 
tlirouchout  the  (  y)  plane  following  the  procedure  described  in  Section  ll.b 
The  sim.ilaritv  par.imeter  values  '  (wb.ich  give  equal  i  nc  rer.en  t  .s  of  t  ai-.g  en.  t  i 1 
veloeitv  for  a  specified  number  of  boundary  layer  points')  .ire  ca  1  e  u  1  a  t  ed.  fro; 
Equation  IS  using  'lewt  on-Raphson  iti-ration.  The  lani  i  \'a  1  on  t  wilm  -  .m  ;; 
given  >;  (  line)  ari.  found  from  the  similaritv  relation,  Ilyuatio::  la. 

derivativt'S  gi\’en  bv  llqinitlons  22  and  2'3  are  then  appro.'-:  ima ;  e.l  ii.sin.  the 
above  value.',  for  >•  and  tlie  specified  ’  diffiTences  a.s  follows: 


■43 


! 


f 
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-  ut  •  =  n 

( -  -  ) 

(jt)'  =  Cy,''  ■  at  •  =  -  1,': 

where  =  1,  y.  are  the  boundarv  layer  y  valuer,  and  'v  is  tiu  syieitiLd 
tar  field  y  difference  between  the  last  tw'o  points  on  the  eiven  '  line. 

Then  tlu-  dat’.pi[ig  factor  D  is  found  using  Kquation  2s.  hext  tlie  ana '  t  ;  c;;  1 
expression  for  d’  /dv  in  Equation  21  is  evaluated  on  the  far  field  inter’cal 
2  \-  and  on  eacii  boundarv  laver  2a-  interval  deterr.ined  hv  ti.e  .  \  alia.  . 

...  .  -  j 

Ttiis  systor.  of  equations  for  (  O  is  solved  using  Cdiuss  el  in.inat  ion. .  I:. us, 

for  a  gi\’en.  '  line,  values  for  and  P  a.re  calculate-d  for  i;sa  ii'.  i  \’a ue  t  i;- 


the  Q  (  g  ’  )  f  ore  in 

g  function.  Tu 

is 

f  unc 

t  ion  is  then  used  durin,-  th 

C'  S'" 

iteration  for  the- 

transf  ormat ion 

val 

UcS 

of  X  and  y. 

Additional  ' 

contours  can.  be 

ad 

d  t.  d 

in  tile  wake  region  by  using 

1  !  K 

interpolating  technique  given  in  Appendix  C. 

B.  N_a_v_ier-Ptoke.'-  I'initL  hifference  iscuations 

Approxiitat >-•  solutions  of  the  unsteady,  Reynolds  averaged,  ir.cor.’prts- 
sible,  turbulent,  Na  v  ier-S  t  okes  equations  are  obtained  by  using  a;- 
finite  d  il'ferenci.'  procedure.  Explicit  nurnerical  tnethods  car,  ;  v  > 

obtain  approxin.att'  solutions  by  advancing  the  tine  by  sn.all  increicA.i;  . 

The  allowable  tiiv.e  increments  are  constrained  by  numerical  sta'  ilit;  ciaa.  rh:, 
lr:;licit  finite  difference  methods  can  also  provide  approx  im.a  t  >.  so’.-;:;,: 
inarching  in  tint,  wit'n  small  time  steps.  Implicit  moth 'dr,  usual'-  d,  ;  ... 

numerical  stahilitv  constraints  imposed  on  the  time  step.  hcwi  .  .  i  ,  ,t:  ...  - 
time  step  a  lar,-,<  sv.-;tem  of  simultaneous  difference  ecpia.tiorr-  ru.-t 
;\  r.atrix  itirat  iv.,  technique  s-uch.  as  SOK  iteration  is  ar.  atlra,  :  iv..  ’.a;  i  - 
rnethod  hi  .  ause  of  the  larg.  m.atrix  sime.  Convergcnci  iri;.  ri.i  .ir.  l'  ;  i-.'i  - 

due  ed  in  place  o!  stal'ilit'.  criteria. 


Ihc'  inpHcil  "lelhod  consists  of  hot!’,  a  linearization  of  Ihe  cqnat  io;.: 
and  the  application  of  finite  differences  to  those  equations.  'ihe  x  co-.- 
ponent  ,  trails  fvirr.ed ,  rao.iientuzi  equation  given  by  Kquation  48  is  con.sideri.  d 
first.  Rc-ar  ran.- inc  the  convective  terrt.s,  the  equation  beci'Cics 
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.'.'quaticn  Sa  is  nov  vriltcii  in  (juas  i  - 1  inea  r  for::;  vlu.  ri  lover  order 
der  i'.-at  ives  of  the  variables  occur  in  "coe  f  i  ic  ienf  s"  .  dhe  current  Velxe 
feir  tliC  unkne'vn  varial-K-  at  location  (i,  j)  occurs  in  the-  iiighest  order 
derivative'  of  e-ao'i  t  -rr.  vi.ile  th.e  ’'e'oef  f  ic  lent  s"  arc  evaluated  wit':',  t'ae-  lo..-t 
available  v.ilue.s  for  tlie  variables  u.  v,  and  p.  In  th.is  forra,  tliC  coef  f  i  ,■  i  eot  t 
of  u  .  in  tlie  convecti'Ce-  terms,  (uy, -vx,  V/,1 ,  is  proportional  te'  lIk-  conpe'nent 
of  ve-locity  in  tiu-  ’  direction.  .Similarly,  the  coefficient  of  u.  in  the  ce'n'.'ec- 
t  i'.'e'  terr.s,  (\'x  .-u;.  d /,i ,  is  p rujuir t  lona  1  to  the-  conponent  of  velocit-.'  in  tin 
'  dire'Ctie'n.  Tiie.se  quantities  are  designated  !'('  and  VC,  respec  t  ive  1  v  . 

Finite  diffe-rence  exp res.s ions  ..-iven  in  .-‘ippendix  B  are  used  to  appix'xi- 
n.itt  t!;e  ter,";.'  in  Kquation  be,  first  order  accurate  hackvard  diff.-ren.e- 
are  u.se-d  feT  tlie'  title  derivative.  Tile  fir.st  derivatives  of  pressure-  are 
app  r  'X  in.itee!  witij  .see'ond  order  accurate  ce-'ntral  d  i  f  f  erence-s  .  1  i  1  e-v  i  se  , 

e  .  nd  einie  r  central  difference'’  express  i  e'n.s  are  u.'od  for  the  .second  de  i"  i  \-,i  t  i -.e  ; 
s!  \-e]e>eit;.-  in  tiie  viscous  terms.  .''e-cond  order  accurate-  upvin.el  d  i  f  !  e  re  n,- i  n  ■ 
i-  u  e-'i  for  tiie-  vele'eit-.  fir.st  derivative.s  in  the-  ceui\'e-e- 1  i  ve  t  e  rr 
qua-'titie  "i  ,in,!  V(  ele-te  rr.i  lU'  tiie  .'iqn  e'f  t  lie  leval  veloiity  in  ti.e  '  an.’  ' 
diroct  ic;n-  .  It  'T  t  lu.-n  hacki-.'ard  iipvind  d  i  f  f  ere-iu' i  n  c  is  usi  si  u-.;-  tin 

.!e  r  i  . ,  1 :  i -.e  v'  ill-  'T  '■  (iirtates  forvarj  upwind  d  i  f  f  e-re  ne  i  no  .  .-M.-'v  .  i.ek- 

A-trd  ip'.-in!  d  i  ‘  ;  er. -a.  i  no  for  '  de-r  i  va  t  i  ve--  oCeUirs  if  \T  !'  w’  i  1  i 


upviniJ  d  i  f  f  <.'rt;nc  inp  is  used  if  VC  '  0.  Seco;u!  onUr  .iccurait  up-iM:  d  i 

ferenees  ar>,  also  used  for  tlie  velocity  first  deriv.ilives  ii:  viscev..- 

terir.s.  The  direction  is  deterr.ini  d  hy  tiic  sign  of  the  coefficient.'  . 
Points  on  the  ,I  =  2  and  I  =  JMAX  -  1  contours  use  fir.-t  ordei'  upwir.f. 
differences  for  ’  first  derivatives  if  tl.e  velocity  dire.tio!:  ir.ditati.' 
flo'.c  fro:-;  the  body  ..iurface-  (J=l')  or  the  outer  boundar-.-  i' J  =  i  .  henr 

the  body  surface  th.e  grid  spacing  is  extreniely  str.all,  and  velocite  cradie 
are  small  near  the  outer  boundarv.  Thus,  in  both:  cases  the  ass.sialed 
artiiicinl  viscosit;.'  has  a  negligible  effect.  The  t  r.-n.-format  ion  deri'.'e- 
tives  in  hquatiot:  89  ara  evaluated  with  second  order  central  d  i  f  f  e  r  ls.c  ■  . 

Tlie  velocity  derivatives  in  the  velocity  divergence  tern,  are  also  e'.'aliiat 
using  central  differences. 

Piiferenco  eq'iiations  are  written  wit'n  the  followinc  con.Ven  lent  co:.- 
ventions.  Terms  which,  contain  unknouT.s  at  location  (i,j)  are  f-allv  s-d-- 
scripted.  tU'  finite  difference  expression  for  location  (i.  '1  a.s 

descri!''ej  previously  is  n.nderstood  when  only  the  torn,  is  .tiven..  The  su:  - 
.ccript.-.  (i,  j)  and  super  .''cr  ipt  n  arc  assuned  when  omitted.  hct.-  result';:-.,- 
ci;:er(,'nce  equation  for  Equation  Sq  be>cor.<.:s 
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for  VY  > 
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-  1 

JVl 

=  j  +  1 

JV2  =  j 

_  2 

JV2 

=  j  +  2 

The  difference  equation  for  conservation  of  mass.  Equation  51,  for  location 
( i ,  j )  is  given  by 


°  =  2J 


u  . 


i-l) 


2J  ^'''i+1  '■’i-l^  2J  ^^j  +  1 


2J  j  +  1  j-1 


(^1) 


The  following  finite  difference  definitions  are  introduced  for  the 
transformation  derivatives  and  coefficients  at  location  (i,  j): 


XF.TA  =  X  /2J 

r, 

x>;i  = 

X./2J 

YETA  =  y  /2J 
r, 

YXl  = 

y,/2J 

1 

7 

T 

ALFA  =  'JJ" 

GAM.\  =  y/J*' 

BETA  = 

-?/2J“ 

With  these  definitions,  Equation  90  can  he  written  to  provide  an  expression 
for  u? .  given  by 


47 


u..  +  lll-YETAU>.^^  -  p._j)  +  V>:i(p.^,  -  p,_^) 


+  7; - 1A1VA(u._^,  +  u.  ,)  +  GA:-PA(u  +  u.  ,) 

Re^  1+]  1-1  j+1  j-1 


+  Bi''rA(u,  -,  ,  -  u .  ,  .  ,  +  II .  ,  .  ,  -  u .  ,  .  , )  +  rv  (Au,,.,  -  u,,,_  ) 

i+1  j+1  1+1  J-1  1-1  J-1  1-1  J+1  1\:  I . 1 

+  vv  /  1  +  rc,  +  VC  +  rv  +  vv  ) 


+  7; -  (ALFA  +  GAMA,')  ] 
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wL.ere  *  dc-notes  tht'  value  of  the  unknoum  calculated  directly  fro-  ti'.c  dif¬ 
ference  equation. 

Several  observations  are  rr.ade  concerning  tht-  difference  equation  for 

given  by  Equation  92.  The  unknowns  v_  and  do  not  appear  explicitly 

In  fact  does  not  appear  at  all.  Tlie  y  component  of  velocity  .  does 

appear  implicitly  in  the  qiiasi-linear  coefficients  IC  and  VC  as  does  u... 

^  J 

The  previous  value  for  u^.  occurs  explicitly  in  the  term  containing  the 

velocity  divergence  D.  The  finite  difference  linearisation  assumption 

requires  that  the  latest  prior  values  for  u..  and  v..  are  used  in  these 

ij  ij 

instances . 

In  an  analogous  manner,  the  y  component,  transformed  momentum  equation , 

Equation  is  written  for  .  and  becomes 

ij 


V ,  .  =  v ,  .  +  At  [  >:E1A(p  -  p  .  )  -  XXI  (p  -  p  .  ) 

11  11  '1+1  '  1-1  '1  +  1  '  I-l 


+  VC  (Av^^,  -  v^^^,')  +  AT  (Av^^^  -  v,^,J  -  V!1 
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+  lALiA(v,^^  +  v.^^)  +  GAMACv  +  V  1 
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1+1  j+1  1+1  j-1 
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Again,  does  not  appear  either  explicitly  or  implicitly.  Trie  unknov,-n 

appears  implicitly  in  the  coefficients  I'C  and  VC.  Also,  occurs 

implicitly  in  I'C  and  \'C  and  explicitly  in  the  velocity  divergence  term. 

Ihe  rriost  current  previous  values  for  u,.  and  v..  are  used  in  these  lower 

ij  ij 

order  coefficients. 

The  static  pressure  throug'nout  the  flow  field  is  calculated  f  ror:,  the 
Poisson  pressure  equation.  Equation  50,  derived  in  Section  II.  Second 
order  accurate  central  differences  are  used  to  approximate  the  .second 
derivatives  of  pressure  in  the  transformed  Laplacian  term.  Tlie  fir.st 
derivatives  of  pressure  in  the  Laplacian  are  approximated  by  second  order 
accurate  upwind  differences.  The  direction  is  based  on  the  sign  of  t)::- 
coefficients  c  and  t  which  is  the  same  method  used  for  the  r.omentum.  equj- 
tions.  In  this  equation,  a  nonlinear  source  term  occurs  which  i.s  comj'osc-d 
of  transformation  and  velocity  derivatives.  The  transformation  der ivat ivws 
are  approximated  with  second  order  accurate  central  differences.  Sccor.d 
order  accurate  upwind  differences  are  used  for  the  velocity  derivativL s . 

The  direction  of  the  differences  at  each  location  is  determined  by  tin- 
direction  of  the  local  ’  and  ’  velocity  components  related  to  VC  and  VC 
previously  defined.  The  time  derivative  of  the  velocity  divergencL  is 
approximated  with  a  first  order  accurate  backward  difference.  iht.  require¬ 
ment  that  the  velocity  divergence  at  the  nth  time  ste;  be  zero  is  incor¬ 
porated  into  the  difference  equation  by  setting  to  xt-ro.  "I  lu  nurorical 


I 
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non-zero  quantity  ^  is  retained  as  a  correction  factor.  In  this  v.-a'.-, 
the  static  pressure  is  found  which  tends  to  satisfy  mass  conservation  at 
every  location  for  each  time  step. 

The  following  definitions  are  given  for  the  finite  difference  repre¬ 
sentation  of  the  source  terms  as  described  above: 

rX  =  (>vu  .  -  y  ,Ai^  )/J 

L'Y  =  (.X 

VX  =  ( V  V  .  -  V  ,v  1  /.] 

■  r  f  ■  '  f; 

\’Y  =  (x  -  X_  v  .)  /.] 


If  th.ese  definitions,  the  previous  transformation  coefficient  definitions, 
and  index  conventions  are  used,  the  transformed  pressure  equation,  Equa¬ 
tion  V',  i^  rearranged  for  the  field  static  pressure  p^ .  given  by 


+  +  2(V.X)(rY)  +  (VY)“  +  AI.F.4  +  d.  ,) 

I  1+  J  1 - i 


^  Pj-1^  BETA(p.^j  -  p.^j 

:-!  ^  ■  PlV2^  ^^Pjvi  ■  P,1V2^  ^ 


^  2  i  ALFA  +  G/V'l/\) 


In  thi'..  equitien  p_  does  not  appear  implicitly.  The  velocity  con',pone::th 


u.j  an!  v_  are  onl;.  present  implicitly  in  the  nonlinear  source  terr.h 


Tte  tliree  difference  equations.  Equations  92  thru  91,  form  a  set  of 
equations  for  the  flow  variables  u,  v,  and  p  at  location  (i,  j  1  in  ti.u'n’.s  of 
the  transformation  derivatives  and  values  of  the  variables  at  neigh! orinr 
points.  Tfie  quasi-linearization  and  differencing  metliods  have  (iiooupltd  tli 

ct  to  the  higher  order  terms  for  the  unknowns  at  loc.ui 


90 


equations  witli 

(  i  ,  j  )  • 


where  s  is  the  iterate  number,  *  designates  the  current  value  obtained 

from  the  difference  equation,  and  c  and  are  acceleration  parameters. 

^  uv  p 

The  acceleration  parameters  are  always  unity  for  the  first  iteration. 

The  computation  of  the  wall  static  pressure  is  carried  out  separately 
as  discussed  in  Section  II.  The  pressure  iteration  technique  of  Chorin  (bS 
given  by  Equation  40,  avoids  the  difficulty  of  formulating  one  sided  dif¬ 
ferences  at  the  body  for  the  Laplacian  of  pressure.  Also,  at  the  body  sur¬ 
face,  conservation  of  mass  is  directly  imposed.  Second  order  accurate- 
upwind  differences  are  used  for  the  transformation  and  velocity  first 
derivatives  in  the  velocity  divergence  of  Equation  51.  T!ie  acceleration 
parameter  1,  which  is  related  to  SOP.  iteration  of  the  Poisson  pressure  c-tpia 
tion,  is  given  by  Hodge  (53)  as 

'  =  2:..  ,  /  .At(ALFA  +  GA.VA)’ 
pb 

where  is  an  acceleration  parameter  for  the  corresponding  StH.  iteration 

of  the  Poisson  pressure  equation.  For  consistency,  the  iteration  e'.aa. i  on 
ir  repeated  a.s 

(s+1)  =  -  :i)  ('-71 

where  s  represents  the  iterate  number  and  the  indt  not.ilio:^  is  ke;)t  as 
described  above. 

The  trailing  edge  pressure  on  an  airfoil  can  he  diflioilt  to-  obtain 
using  the  iteration  procedure  if  a  sufficientlv  fine  yrii!  is  vu't  iires.nt. 


Hodge  (53)  experienced  such  difficulties  in  his  lair-inar  flow  work.  The 
trailing  edge  pressure  can  then  be  calculated  using  an  average  oi  value- 
obtained  by  linearly  extrapolating  the  pressure;  calculated  at  the  closest 
points  on  the  upper  and  lower  surfaces.  For  a  two  point  linear  extrapola¬ 
tion  in  the  y.  direction,  the  pressure  at  the  trailing  edge  becor.LS 


(X]  -  x,^)  p^  -  (x^  -  p^ 

2(x^  -  x,-^) 

'>(•.•  -  X  i 


( ^  ~  a  ) 


where  the  '  s  are  the  x  values;  at  tht  selected  surface  points  an.d  tlie  sur¬ 
face  j  =  1  value  is  understood. 

The  iteration  procedure  follows  a  prescribed  order  at  each  ;K.'int  loca¬ 
tion  (i,  j).  Only  the  body  surface  pr'^ssure  is  cortputel  on  hod;.'  points. 

At  field  points,  the  u  velocity  component  is  coTtputt'd  first,  ti.en  t'ne  v 

velocity  component,  and  finally  the  field  pve,ss\jre  p.  Ihc  point.*-  uri.  ..  rdi  i'l  . 

by  varying  the  '"  ( j )  index  from  j  =  1  thru  j  -  .lMAX-1  or-,  succes  s  i  vi  ’  litee 
in  the  tran.sformed  plane.  The  ’  ines  i  =  1  thru  i  =  IM-Al  -]  ari  t  ra'.'i.  r  ed  . 

Values  for  the  variables  at  i  =  1  are  set  equal  to  tht  'calues  at  i  =*  I 

to  insure  continuity  across  the  t  ran.-^  f  orma  t  ion  bran;':,  cut.  The  i  ir.di 
is  adjusted  across  the  branch  cut  in  a  manner  identic.;]  with  th.ut  dose  r  i  !-.i  <: 
for  obtaining  the  transformation  fr'ir  equations  FT  a;'.;’ 


C.  Bonn  ' .ary  and  Initial  (hyncH tj on -- 

Boundary  and  initial  conditions  are  formnlatid  for  usi  in  the  eonpuia- 
tional  plane  in  order  to  completely  define  I  lie  numeriral  prohle:  .  hiuKular-. 
conditions  on  the  hod;,  surface  are  conven  ieiit  1  v  specified  with  tin  li  '  - 
fitted  c  c'ord  i  na  t  e  svstem.  The  first  '  contour  is  con-truite.!  te  '  <  the 


i 


body  surface.  The  analytical  no  slip  boundary  conditions  for  velocity 
given  by  Equation  52  can  he  directly  applied  to  the  body  surface  in  the 
computational  plane  and  are  given  by 

u  =  0  V.  =  0  (98) 

b  b 

where  the  subscript  b  denotes  the  body  surface  contour.  Tlie  surface 
pressure  is  calculated  as  discussed  in  Section  IV. B. 

The  lines  i  =  1  and  i  =  IMAX  in  the  computational  plane  define  tin- 
branch  cut  in  the  physical  plane.  Ko  boundary  conditions  are  required  c^r 
allowed  on  these  contours.  The  matching  of  solutions  on  these  lines,  wl.icl 
was  described  previously,  provides  a  cont  intiou.s  .solutio.n  a.-ross  t’le  cut. 

The  outer  boundary  in  the  physical  plane  frequently  represent.s  con¬ 
ditions  at  infinity  (freestream  conditions)  given  by  Kquaticn.s  53  and  54. 

In  the  body-fitted  coordinate  system,  the  outer  lioundarv  becomes  the  largest 
line  in  the  com.putat  ional  plane.  T'ne  physical  o\iter  ’'oundar'.-  i-  those:; 
to  be  a  circle  with  a  radius  of  10  airfoil  chord.s.  Die  d  i  s  t  r  ibu  t  i  cm.  of 
points  on  the  outer  boundary  is  readily  specified  t'n  a  circK;  al-o,  il., 
initial  input  solution  for  obtaining  the  transfort.a! ion  a-  discus-ed  in 
Section  IV. A  is  simple  to  specify.  Ghia  and  Hodge  (8.1'  applii  1  a  r  i  o.,  ’ 

inviscid  anlaysis  to  a  .ioukowski  .airfoil  at  angle  of  .Uiack  usin,-  fri,<- 
stream  boundary  conditions  fc-r  different  outer  boundarie.-- .  Tliey  foin'd 
tt'.ut  for  a  10°  angle  of  attack  the  lift  coef  f  ic  a  et'.t  differed  .  Dst  t’l.n; 
one  percent  and  the  maximumi  suction  pressure  coe  f  I  ic  it-r.  t  by  .d-out  twe-  [u  r- 
cent  frov,,  the  analytical  result.  ha  vie  r-.*^  t  okes  solutions  usiiig  a  circular 
outer  boundary  with  a  smaller  radius  are  compared  wit’  the  It  chord,  r.eliu^ 
result  to  dettrm.Lne  effects  of  houndarv  plucemei-.t. 

With  the  outer  houndarv  elected  to  .-jppr  ox  i  ma  I  i  !'■  far  fii  Id  1  n  .  - 
stream,  tlie  boundary  conditions  for  velm-ity  ;uul  pressure  ari.  spci  i;ied. 

'it 


As  pointed  out  by  Roache  (83),  caution  must  be  exercised  in  applying  tiie 
analytical  conditions.  Equations  53  and  54,  which  are  strictly  valid  only 
in  the  limit  of  large  distances  from  the  body.  If  these  conditions  tre 


used  on  a  boundary  at  a  finite  distance  away,  the  result  may  predict  no 
drag  since  a  wake  cannot  exist  at  the  outer  boundary.  This  problerr.  is 
avoided  by  applying  upstream  and  downstream  boundary  conditions  on  the 
outer  boundary.  The  upstream  boundary  along  which  conditions  are  fixed 
is  defined  by  the  semicircle  in  the  half  plane  x  0  where  the  x  axis  lies 
along  the  airfoil  chord  with  the  origin  at  the  midchord.  The  downstrea::. 
boundary  for  which  some  variables  are  permitted  to  be  free  is  the  semicircle 
in  the  half  plane  x  >  0. 

The  upstream  boundary  conditions  for  the  incoming  undisturbed  flov; 
become 


u,  =  cos  u 

1  J'lAl. 


V,  =  sin  j 

1  JM\.\ 


»!  j'la;.;  ■  “  «« 

where  i  ranges  over  all  T  contours  which  terminate  on  the  upstream,  boundary. 

The  downstream  boundary  conditioT s  must  allow  a  wake  dowmistream  of  a 
body  to  pass  through  tVie  boundary.  The  no  change  boundary  condition  has 
been  widely  used  (53,83)  and  meets  this  requirement.  One  approacli  requires 
no  change  of  the  velocity  components  in  the  freestream.  direction.  I'.'-ing 
the  expressions  for  the  directional  derivative  and  gradient  in  the  trcins- 
formed  plane,  this  boundary  condition  can  be  written  for  the  u  component 
of  velocity  as 

u^  =  u (x^  sin>  -  y  cos*)/(x.  sin'  -  cosO  (I'i,', 


where  o  is  the  geometric  angle  of  attack.  Mow  ron.sider  a  .similar  chnen- 
stream  boundary  condition  where  the  no  change  condition  is  in.po.- eii  on  tla 
velocity  components  along  the  down.stroam  ’  contours.  In  a  sirwlar  mannir, 
this  condition  in  the  transformed  plant  becomes 

ip  =0  v^  -  n  (I".' 


This  boundary  condition  is  identical  with  the  freestream  direction  bounde.ry 
condition  of  Equation  100a  when  the  y  contours  are  in  the  freestrear.  direction. 
In  this  case,  tan u  =  and  Equation  100a  becones  identical  with  Equation.  1'' 

This  case  does  occur  in  the  downstream  wake  region  for  the  transformation 
obtained  in  Section  IV. A.  Small  changes  in  the  velocity  components  are 
found  in  the  far  wake  region  of  bodies.  The  downstream  far  field  invisciu 
flow  field  also  has  negligible  changes  as  seen  in  the  analysis  of  .Aj.pend i.x 
r.  Thus,  the  no  change  conditions  of  Equation  100b  approximate  the  sm.all  far 
field  dow'nstream  variations  in  velocity  while  allowing  the  momentum  defect 
in  a  wake  to  exit  at  the  outer  boundary. 

The  downstream  outer  boundary  condition  for  pressure  must  he  considt^red 
along  with  the  velocity  conditions.  Hodge  (53)  investigated  the  use  a 
similar  no  change  condition  for  pressure,  namely,  p^.  =  0  coupled  witli  the 
velocity  conditions  of  Equation  100b.  His  laminar  solution  developed  pres¬ 
sure  oscillations.  The  more  restrictive  condition  of  specifying  the  free¬ 
strear.  pressure  in  the  far  wake  was  successfully  used  by  Hodge  (53').  th.e 
specification  of  pressure  recovery  witii  a  freely  developing  velocity  defect 
in  the  far  wake  gives  a  set  of  numerical  boundary  conditions  wViich 
can  physcially  represent  the  far  wake  at  a  finite  far  field  dowist  rea::; 
boundary.  The  downstream  pressure  boundarv  condition  becomes 

The  suitable  description  of  the  physics  at  Itie  downstreati  comjs'.ta- 
tional  boundary  must  be  written  in  finite  difference  form.  The  no  cii.is.m 
condition  approximated  with  a  second  ordc-r  accurate  central  d  i  l 1 1 ren , a 
applied  at  the  midpoint  in  the  computational  plane  lietweon  tho  r- 

point  and  first  interior  point  on  a  downstream.  '  contc-ur  (i  index'  ci\i 


'‘.ELI--:  =  ‘'.iiEd.:-] 
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for  the  downstream  boundary  conditions  of  the  velocity  cor.;ponent  s .  .'-.r. 

alternate  forinulation  Is  obtained  by  approximating  the  no  change  condition 
with  a  second  order  accurate  upwind  difference  applied  at  the  boundary  point 
Then,  the  douTistrean;  boundary  conditions  for  the  velocity  corgioiu  nt  s  hecon.L 


'■'jMAX  ^^'■JM\X-1  '■'j'TAX-dl/j 


These  boundary  conditions  also  follow  from  a  second  degree  polynor.inal  fit 
with  a  zero  derivative  imposed  at  the  downstream  boundary. 

The  time  dependent  finite  difference  method  given  previously  recniri' 
a  Set  of  initial  values  for  tlje  velocities  and  pressures  at  all  tin,  H.  ; ) 
locations.  This  requirement  parallels  the  result  discussed  for  the  dif¬ 
ferential  equations  in  Section  lll.b.  Two  methods  ar<.’  used  to  provide 
initial  conditions.  In  the  first  method,  an  invisoic!  solution  will,  cir¬ 
culation  for  the  flow  field  at  the  given  angle  »'f  attack  is  obtained  usin,:.', 
SOK  iteration  as  described  in  .‘.ppendix  G.  This  solutioii  then  pri-\'ide,- 
initial  values  for  the  havier-Stoke.-.  solution.  Tiie  second  techniqvK  usls 
the  Navier-Stokes  solution  for  one  angle  of  attack  and  rotute.s  thie  n’i  1  c-c  i  i  i  i 
by  the  ch;ing,e  in  ancle  of  attack.  TIk-  si't  of  in.itial  valuis  for  vlU  ,.  if  on 
pressure  is  then  used  to  start  a  Xav  ier-S"  oki-s  solution  at  tin.  new  an.i-lu  of 
a  1 1  a  c  k  . 


li .  Tin'huJj.'nce  _Xyydtyl 

The  numerical  procedure  wliich  implenienls  the  eddy  visco.-,itv  turhulxn.ci 
model  described  in  .Section  III  i  ,s  discu.ssed.  Ti.e  two  la.'er  mod.  el  r.e.ii  the 
airfoil  surface  is  presented  followed  hv  t!ie  far  wake  modi  1.  ihe  pinm.dnri 
used  in  the  near  w.ike  is  then  described.  The  edib.  viscosit-  di^;  rihn!  i  'n 
throughout  the  flow  fit-id  is  taleulated  at  t  lu  beginning  of  eaf'  t  i:  m  !  -  ;  . 


SOR  iteration  is  then  used  to  determine  the  velocities  and  fire.-=surcs  at  thi^ 


new  time  level.  The  eddy  viscosity  distribution  is  not  changed  during  the 
iteration  process. 

The  inner  viscositv  computation  from  Equation  6"  requires  values  of 
the  tangential  and  normal  velocity  derivatives  and  the  distance  fro;:,  th^ 
surface.  These  velocity  derivatives  are  given  by  Equations  70  and  71  in 
the  transformed  plane.  They  are  evaluated  using  second  order  accurate  central 
differences  for  both  the  transformation  derivatives  and  the  u  an.i  v  velocity 
component  derivatives  with  respect  to  '  and  ' .  The  inner  eddy  viscL'sity  is 
identically  zero  at  the  surface  (’=1)  and  need  not  be  calculattd  .sey..r;:'.e’.  y . 

The  corresponding  velocity  expression  which  is  found  in  t'ne  ex  ton  o;'  the 
damping  factor,  Equation  62,  must  be  evaluated  at  the  airfoil  surface  far 
eacli  '  contour  in  regions  of  turbulence.  The  derivatives  with;  resiuct  ti 
^  are  approximated  by  second  order  accurate  forward  differences  found  i 
.-ippendix  B,  and  the  derivatix’es  with  respect  to  ■’  are  evaluated  using  second 
order  accurate  central  differences.  The  distance  fro::;  the  surfai  l  t  •  ti.e 
point  (i,  j)  is  measured  along  the  con.stant  '  line  reprisented  !;.■  i::de:-: 

Tlie  freestrear.  chord  Reynolds  number  appears  throughout  thie  mod<.  1  a 
result  of  w'ritlng  the  equations  in  nondimens  ional  form  and  is  .an  input  par,:- 
met  er . 

The  adverse  pressure  gradient  turbulence  model  .mod i ' icat ion  in  tin 
tr.iiling  edge  region  give:;  by  Equation  (1.0  is  imp]  one:  t  .  ;  next.  '”.1  .ri.  h 
of  the  distance  s,  measured  alone  the-  .airfoil  .'.urface  1  ro: .  thi  ;  r.;  ; i’:  ■  id..', 
separation  point,  is  determined  bv  t’x:;::  i  n  ing  tht-  u  vtlo.  it'.’  co::  po: , .  : .  t  ,1’  p  in: 
on  t’le  ■  contour  next  to  the  urf.;.!  ('“ji.  l!n  se.arih  iegii.  'x:  ]!;■■ 

corresp  :;din.'  to  tiu-  ‘■'0  f)t  rent  lo,  a;  iun  .;nd  ad'.'.'i  cio  tow.iix.  ■li.i  .lir:  :i 

leading  e.l.o  alone  the  sei'oiu:  '  li  ■.  '  ht  ‘.-ar':,  t  er-.  i  ::a  t  .  x,  wh.  :.  .:  jO'  :  ;  i  ,  . 

v.i  1  uc  for  u  .  ,  i  d. ' '  I  i-d  w:  i  '  i  •:  e  it.  a  t  t  .o  Ui  -d  !  1  ow  .  "  he  r v  1  ;  :  i  o ■ ,  •  :  ;  1 
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f  is  then  evaluated  for  eacii  applicable  •  lint-  w!;ere  trailing:  ed^L-  separa¬ 
tion  is  detected.  Tin.-  inner  eddy  viscosity  is  r'.ul t  ipl  i ed  by  the  r e- 1  axa l  it  :: 
factor  which  adjusts  the  adverse  pressure  gradient  inner  eddy  viscosit;.' 
constant 

.\fter  the  inner  eddy  viscosity  is  computed,  the  limiting  tecVinicue 
evalutiten  the  computed  value.  Tlie  inner  eddy  viscosity  is  prevented  f ro:r. 
decreasing  in  the  outward  normal  and  downstream  taneential  flo'c  directions 
by  comparing  the  value  with  previously  computed  values  as  gi\'en  by  t'le 
following  sequence: 


ij  -  i  j-1 


jj  -  i-l  j 


T'r.e  larger  vtilue  is  selected  during  each;  comparison..  Ti,e  initial  in::e;' 
eddy  visco.sity  limiting  distribution  is  calculated  in  the  attacbied  hour.dar;.- 
layer  region  with  a  favora'.'le  pr.-ssure  gradient  near  tbie  airfoil  leading  ed. 
The  first  or  outward  direction  limiting  condition  is  impo.-.ed  here  .:i.s  well. 

The  limit  value  of  the  inner  eddy  viscosity  is  then  compared  wiibi  tin. 
calculated  outer  layer  edd;.'  visco.-itv.  V.'h.enc-ver  tV;e  inner  edd;.'  vi  .sco.s  i  t 
first  e.'vceed.s  tlie  outer  value,  the  outer  eddy  viscosity  i  .s  then,  used,  fcr  tin 
rer.'.aining  locations  O'utward  from  the  surface.  Tliis  switch  insr.res  Ih^  con.- 
tinuity  of  tb.e  eddy  viscosity  d  istr  iSut  io:'. . 

The  computation  of  the  outer  eddy  visco.sity  given  by  Kquation.  T1  re¬ 
quires  values  for  the  local  boundary  ]a;.er  edge  velcuity  and  the  d ;  s;.T  a.  c::.e: 
thickness  defined  by  brpiation  64.  "iise-se  quantities  arc  calculated  for  i  ’ 
contour  whiobi  crosses  the  turliulent  boundar-,'  layer.  Tor  a  '  line,  tin  Km,.' 
edge  velocity  u  is  defined  to  be  the  mat-:  im.ur.  tancentiai  '.'l- 1  o:  i ! '■  r.i  i.  uved. 


rclativi  to  tlie  local  siirfact-  tangent  line.  'Wv  tangent  i.il 


>  '  1  e  ■  ■  i  1-  i'  ;  I 


b;.'  Fqiiation  6b  and  jy  evaluated  as  describe:!  pre'yiously  tor  t’u  i' on  r  ed,.' 
viscosity.  The  houndar'.  la'.'er  thickness  '  is  llieii  tbi  distance  .ro'i,  lb. 


lir.  fret;  the  surf.t.'i  to  the  first  point  win  re  tin  tanginti.'i!  y;  bo  it 


Ifss  than  u  .  'Ihu  displacer-.eiiL  thickness  is  calculated  usiae  traja- 

e 

zoidal  integration  where  the  liit.its  of  integration  art'  de t err.i Tu.-d  ti  e 
computed  boundary  layer  thickness. 

Tile  resulting  eddv  viscosity  distribution  is  neat  modified  tc  simulate 
the  transition  from  laminar  to  fully  turbulent  flow.  iae  transition  factor 
r,  given  by  Equation  67,  is  calculated  as  follows.  'i'h.e  startinc  loratii.n 
for  transition  on  th.e  airfoil  surface  must  be  specified.  This  locatioi:  is 
specified  by  designating  the  first  grid  point  on  the  surface  downstrecc  as 
Well  as  the  distance  to  that  first  grid  point  .',s.  The  transition  di>tar.tc 
s  from,  th.e  starting  location  to  a  given.  '  line  is  calculated  hy  adding 
to  tile  distance  between  the  successive  surface  grid  p^int.s  with  known  loia- 
tion.s.  The  relaxation  factor  '.  is  determined  from  the  transiticr:  factor 
Equation  67  and  the  definition  of  '. .  If  i  =  3/d  is  specified  wlien.  ti.^ 
transition  distance  equals  three  fourths  of  the  total  transition  lengtli,  tlu 
1/.’.  is  given  hy 

!/■  =  2.Hd56/x.  (It.hi 

wliere  Xj  is  tlie  total  assum.ed  transition  length.  Th.e  eddy  viscosit;.  di--- 
trihution  computed  previously  on  a  '  line  is  multiplied  by  tlu  corres;  on,!  in.c 
constant  "  to  obtain  the  revised  distribution  of  lur'.iuler.t  c-dd;.-  ’.■i.-c,-  it  ’ 
wiiich.  Includes  tlic  transition  region. 

The  final  modification  made  to  the  eddy  visco.sitv  d  i  s  t  r  i  h'a  t  i  on,  aero 
the  upper  surface  turbulent  region  involve.?  the  decreasL  of  1 1;  r:  ; ;  ’ .  ■ . .  in 

till  diri'Ction  from  the  boundary  layer  toward  the  far  field.  'i'u  i  n :  >  n:  i  t  1 1  n 
factor  giv.’n,  bv  !!;;uatiou  63  models  this  beh.avior.  'il.i  ]'rev  i  ou:- ■  oa  1 .  ■.; !  ,i  t  . 
hound, iry  la;.er  t !,  i  <  kne.o.o  '  and  tin-  di.staniu  from,  thi  airfoil  surl,.,.  t.  t’u 
(i,  j)  lo.atior.  or.  a  liiu  are  used  to  compute  the  in  I  erm  i  t  ti-nc  v  to;  t:,.w. 
looati  'i-;.  ’1  he  previous  v.ilues  of  the  outer  edd\  viscositv  are  mult  ip!  ied;  ' 

th.  c r  r  I  ■  -  ponl  i  n v.ilui  of  the  i  nt  i-rn.  i  1 1  ency  farti.r  to.  ol'tain  the  iina!  o-,!  h 
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viscosity  d  i  s  t  r  ibu  t  io;-;  for  tiio  turbvileiit  rc-;/iop.  oliovc-  tin.  iirlcii  sur.'oii. 

,  til 

at  the  n  tine  step. 

rile  procedure  for  conputini;  the  eddy  viscosity  distribution  on  n. 
line  acrofls  the  t u rlni  1 en t  flow  zone  above  the  airfoil  is  repeated  down- 
strear.:  lor  each  ’  line  through  the  95  percent  chord  point.  li.e  '  contour- 
near  tile  trailinc  edce  deviate  Iron  the  local  nornal  direction  in  exce.-.- 
o;  1  ih  .  For  locations  in  this  snail  region,  the  eddy  viscosity  is  found 
by  usinc,  the  last  calculated  value  upstrean  in  a  direction  tangent  tr  I’n 
.iirfoil  surface-. 

li'.e  calculation  of  tlie  far  wukc-  eddy  viscosit,.'  given  by  h yu a i  i 7.:  is 
sir-ilai'  to  t:u-  coinput  at  ion  of  the  outer  layer  eddv  visco.sity.  In  tin- 
recic'ii  the  •  contour.s  are  approxir.ately  orthogonal  to  the  flc-'v  directic-n  in 
the  vakc-.  Till.  ’  contour.s  are  thus  conveniently  defined  patiis  acres.-,  the  vaf.c 
for  use  in  co:;. ptiting  wake  characteristics. 

1 'n  an  ’  cont-.'ur,  t!ie  component  of  velocity  in  the  treesZTeon  direct  i.-n 
is  computed,  for  the  grid  points  spanning  the  wake.  Tiie  min  imurr.  of  tl.e.sc 
directional  velocity  components  in  the  wake  is  determined.  Tiic  riac:;:--.: 
v,!luo.-=  of  the  velocitv  components  which,  occur  on  cither  side  of  the  tirim".;- 
value  are  found  and  designated  the  upper  and  lower  e-dge  ve  1  or  i  t  i  e--- .  i'vn 
the  uppt-r  and  Inwc-r  edg;es  of  the  wake  arc  defined  as  tlie  points  wb.erc  th. 
Ic'cal  Velocity  component  in  the  fret-stream  direction  ‘‘irst  reaches  p-..  r- 
cei’t  o-'  the  upper  and  lower  edge  Vt-locitios,  I'c'-spcc  t  i  ve  1  v .  ’he  dist,o-;ci-  •  r, 
the-  y.jint  r<:'  minlmun  velocity  to  tlie  uppe-r  and  lo-we-r  wai-u  c‘d.-i.  ^  .’.ri  ihfi--,-,  d. 

as  tile  ujipe  r  and  lower  wake  tiiickncsses  ,  re.spec  t  i  ve  1 .  These  dist.us  <_ 
d  e  t  err,  i  lU'  tlie-  lir.its  of  integration  when  ca  1  ul  a  t  i  n tbie  upgi  r  .me  Imcer  w.T.n 
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the  average  of  t!ie  upper  and  lover  displaceirient  thicknesse.s .  Ti.e  far 
vake  eddy  viscosity  is  calculated  from  Equation  73  using  the  thickne.s.s  y 
and  an  average  of  t'n-  upper  and  lower  edge  velocities, 

Tt'.e  intenri  1 1  Latov  fttctor  for  each  grid  pc  int  on  an  '  conrtntr  is  deti.:' 
mined  using  the  average  value  of  the  upper  and  lower  wake  t r. i c kne.L- .Se -  f^r 
in  Erpaation  63.  The  local  distance  to  each  point  is  calculated  fro-  the 
location  of  r.iniTnum  velocity  along  the  '  contovsr. 

Th.e  near  wake  model  is  numerically  implemented  by  con.siderinn  tlte  upp 
turbulent  and  lover  laminar  boundary  layers  sepa*'ately.  Th.e  eddy  viscc'.-:i 
distribution,  which  is  calculated  near  the  trailing  edge  of  the  airfoil  .at 
tile  95  percent  cliord  location,  is  extended  into  tiie  near  vake  alor.,', 
direction  parnlltl  to  the  airfoil  surface  defined  by  trie  93  percent  cim’rd 
and  99  percent  chord  grid  point  locations.  The  ciiord  line  is  the  line  of 
extension  for  the  boundary  lavcrrs  in  the  near  wake.  The  ’  lines  i  =  1  and 
i  =  lii.\h-l  define  -he  interaction  zone  in  the  near  vake.  T’ne  non-intL'r- 
acting  di-'tance:  k'L.j  in  Kcuation  74  i.s  based  on  the  lengt'n  of  the  lover 
laminar  boundarv  layer  thickness  computed  at  the  ^‘5  perc  nt  c’.ord  location 
Thi.s  boundary  layer  chickne.ss  is  computed  by  the  same  method  used  p '.-l'-.- i  .".s 
tile  upper  surface  boundary  laver.  The  distance  ’.\T. ,  along  trie  ci.crd 
.  inc  from  tiie  trilling  ed  o  to  th.e  location  vhere  the  far  vaki'  eddy  ■.i.-- 
co^if  modi.!  hegins  is  al.so  calculated  as  a  multiple  of  the  laminar  tr.-.ii- 
inc  .vdee  b  jndary  layer  thickness.  The  eddy  viscosity  in  the  interact  iei' 
zone  is  computed  vith  the  two  calculated  di.  tances  and  tiie  far  wa'-;L  eddy 
viscosit;.'  using  Equation  74. 

E.  Force  (Coefficients 

file  forces  and  moments  exerted  on  the  airfoil  hodv  bv  t  iii  flov,  field 
are  ev.i  1 'ua  t  ed.  at  tin.'  bodv  surf, ace  using  Equation.s  7“,  b'l,  and  b.’ .  i  hi 


integrals  are  evaluated  using  trapezoidal  integration.  The  product  ter:, 
in  the  integrands  of  the  form  (fgl  are  integrated  by  using  product.?  o: 
averages,  (f.  +  f.  ,1  (g.  -  g.  ,)/4.  The  "  derivatives  are  evaluated  u.‘~inc 

second  order  accuratt-  central  differences,  and  the  ’  derivatives  are  forr.;- 
lated  vith  Second  order  accurate  upv;ind  forward  di  f  f  e  re:;ce  ~ .  The  lift  ar..' 
drag  coe  f  f  icietit  .s  are  calculated  by  adding  the  components  of  and  y  directic:. 
force  coefficients  in  the  normal  freestream.  (90  degrees  coun  t  ere  lockwi.se )  and 
f reestream.  directions,  respectively,  using  Equations  Imld  and  I'. 14. 

Th.e  alternate  technique  based  on  the  control  volu::.e  analysis  in 
.Appendi.x  P  for  computing  the  forces  on  the  airfoil  hody  is  num.er  ical y 
implemented  in  a  similar  manner.  The  force  coefficients  in  the  a:id  y 
directions  for  a  two-dimensional  control  volume  defined  by  an  ’  lint-  ai'c 
given  by  Equations  81  and  84,  respectively.  Tliu  line  integrals  art.-  a,c.';i:'. 
computed  with  trapezoidal  integration  using  products  of  averages.  The  ’ 
derivatives  are  evaluated  with,  second  order  accurate  central  differences. 

T!ic  ’  derivatives  are  cvaluaLod  with  second  order  accurate  central  dif¬ 
ferences  on  the  interior  ’  lines.  Seco:id  order  accurate  upwind  d  i  f  re  :'ic  e.s 
are  used  to  approximate  the  ’  derivatives  on  th.e  body  and  outer  bound. ir-.' 

'  lines.  Tb.e  area  integr.al  is  also  evaluated  using  trapezoidal  integration. 
The  Jacobian  for  an  area  element  bounded  by  the  i  and  i-il  line-  .and  ;-l 
and  j  r  lines  is  calculated  by  taking  the  average  of  the  l.acobiar.  v.jln.is 
previou.sly  computed  for  the  viscous  term  in  the  1  inL  integrals  at  tbi  cor¬ 
responding  "  segments  on  the  j-1  and  j  '  lines.  Tim  .irt.i  i:.tegrals  .are 
computed  for  hcjtV.  n  and  n-1  time  intervals  using  tlu'  previou.-l;.  ce-yi.l  ..d 
flow  field  valuts.  The  time  derivative  is  tlien  evaluated  usin.g  ’  \r  l 
order  arciir.ite  backward  d  i ''f  eretn.  e  . 
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Di  scrssio;';  q.-i  rksi'Lts 

N'ur’.erical  solutions  for  two-d ireens ional  incompressible  turb;;]L-nl 
viscous  flow  over  a  NACA  0012  airfoil  section  are  obtained  usinr  t’;e 
implicit  finite  difference  method  previously  described.  llie  '.'ArA  0(>12 
airfoil  section  is  chosen  because  of  its  widespread  use  as  a  test  case 
in  both,  experimental  and  computational  work  (65,66).  Tlie  solutions  are 
for  ancles  of  attack  of  5',  7.5',  9.5’,  and  11.5'  at  a  ch.ord  Reyncld.s 
number  of  170,000.  This  Reynolds  number  is  selected  for  cor.par isoi'.  with 
both  previous  NA(h\  data  (5c)  and  for  data  obtained  as  part  of  tb.is  in'ce.sti 
gation  at  tlie  .Air  Force  F.  d.  Seiler  Research  Laborator;.-  (.'■A,’. 

The  numerically  generated  body-fitted  coordinate  systev,.  w't.ich.  wa.s 
used  for  each  solution  is  presented.  The  ,succe.''S  i  ve-ve  r-rc’ l.ixa  t  i  on 
(SDK)  iteration  numerical  paramet er.s  and  convergimci-  criteria  along  with 
the  steady  state  convergence  criteria  are  then.  discus.<ed.  The  velocity 
and  pres.sure  fields  and  streamline  contour.s  obtained  from,  the  numerical 
solutions  are  examined.  The  calculated  laminar  separation  bubble  charac¬ 
teristics  are  compared  with  empirical  results.  The  position  c'f  tran.-ition 
to  turbulence  and  the  transition  length  relative  to  the  sepesration  bvibVlc 
location  arc^  discussed.  The  computed  airfoil  surface  mean  pre.c.sure  d  i 1  r  i - 
butions  are  compared  with  both  an  inviscid  solution  and  exp.c  r  im.cnt  a  1  data. 
The  num.ericallv  predicted  lift  and  drag  forces  exerted  by  thc^  flow  field  e 
the  airfoil  are  presented  along  with  available  exper  itienta  1  measure  r.i;' .  . 
The  effect^'  t'nat  the  placement  and  type  of  imposed  far  field  boundar-,' 
conditions  have  on  the  numerical  solutimi  are  discusstui.  The  stn.‘  it  ivitc 
of  the  numerical  solution  to  grid  sine  1  .s  examined. 

A  > 


The  curvilinear  body-fitted  grid  shown  in  Figure  3  is  numerical  1;. 
generated  using  boundary  layer  parameters  for  the  attraction  of  '  con¬ 
tours  towards  the  airfoil  surface  and  linear  interpolation  for  wake 
resolution  given  by  Appendix  C.  The  Blasius  series  flat  plate  chord 
Reynolds  number  of  200,000,  a  nondimensionalized  boundary  layer  edgy 
■■  velocitv  of  1.2,  and  seven  grid  points  specified  in  the  boundary  layer 
are  used  in  the  ’  contour  attraction  technique  (7).  The  distribution  o; 
airfoil  surface  points  is  accomplished  with  Equation  26  where  A1  =  0.0, 

A2  =  0.4,  and  A3  =  1.0.  Seventy-one  7  or  surface  points  together  ■■•itl. 
forty-four  '  contours  are  used.  The  transformation  difference  equ.it,  ion-. 
Equations  85  and  86,  were  then  solved  using  SOR  iteration.  Trie  optir-,;::, 
local  acceleration  parameters  varied  between  0.8  and  1.5  throughout  the-  *.  i  uv 
field  with  essentially  no  change  in  the  local  values  after  25  iterati.'rn  . 

After  350  iterations,  the  maximum  errors  in  the  solution  for  botli  and 

-4 

y  were  less  than  ,10  ,  and  occurred  near  the  outer  boundary  above  the 

airfoil.  Additional  7  lines  were  added  in  the  airfoil  val’.e  region  usin,” 
the  interpolation  metliod  described  in  Appendix  C.  F.ight  grid  points 
(four  at  a  time')  were  inserted  on  the  rounded  trailing  edge.  Tabic  fh  1 
gives  the  '  and  ’  indexing  for  the  final  79x44  grid  shown  in  Fivure  1 
and  the  equivalent  indicts  for  the  71x44  original  grid. 

The  airfe'"]  is  thus,  defined  by  a  total  of  79  points.  Eleven  vrid 
points  dt'f  ine  the  first  5  chord  nose  re>;ion  and  17  points  are  used  to  ch- 
fine  the  rounded  trailing’  edve.  Tl-.e  max  imuc,  distance  between,  sur ,  i,- i  ■,«. 
surface  grid  points  is  5  of  the  chord.  The  ’  lines  intersi.'Ot  t'hi  surlaiL 
contour  within  10"  of  the  local  norm.il  direction  exc  ■  j't  in  t  h.e  last  3 
chord  region  near  th.,  trailing  edge. 

The  SOR  iteration  paraneters  .ind  convergence  c  riteria  vchic’i  are  re  .him  d. 
by  the  numerical  implicit  procedure  are  similar  for  eac'n  ancle  att  ni 
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solution.  The  accleration  paraneters  had  values  of  1.0  for  the  velocity 
cor.ponents,  0.9  for  the  surface  pessure,  and  1.10  or  I.IS  for  field  pressure 
The  larger  value  of  was  used  at  9.5'  and  11.5°  angles  oi  attack.  The 
convergence  criteria  for  the  SOR  iteration  required  that  the  maxinuri  change 
in  the  relative  magnitudes  of  u,  v,  and  p  for  all  locations  be  les.s  than 
1,';.  This  criteria  was  relaxed  to  S''  for  the  9.5'  and  11.5'  cases.  Ti'.rec 
iterations  per  time  step  were  the  ma  ;imum  necessary  during  most  of  each 
computed  case.  For  the  9.5'  case,  the  approach  to  the  steady  state  solu¬ 
tion  was  re-run  with  ,  =0.95  and  w  =  1.10  for  two  characteristic  times. 

pb  p 

An  additional  iteration  per  time  step  was  required  to  maintain  the  sac.', 
error  tolerance.  Computed  values  of  lift  and  drag  were  virtually  identical. 

Although  implicit  SOR  methods  have  been  shown  by  linear  stability 
analyses  to  be  unconditionally  stable,  convergence  at  a  given  time  re¬ 
quired  a  snail  time  step  At  0.001.  The  analysis  in  .Appendix  E  based 
on  linear  matrix  theory  indicates  that  the  diagonal  dominance  of  the  set 
of  finite  difference  equations  given  by  Equations  92  anJ  93  has  a  defin¬ 
ite  dependence  on  l/.'.t.  Diagonal  dominance  is  a  necessary  and  sufficient 
condition  for  convergence  (81)  of  constant  coefficient  linear  system.^ 
using  SOR  iteration.  Although  this  system  of  equations  lias  lowc-r  order 
nonlinearities  and  variable  coefficients,  the  order  of  magnitude  study 
of  Appendix  K  indicates  convergence  if  the  most  stringent  requirement 
.'.t  0.CIO03  is  applied.  A  constant  time  step  equal  to  0.0005  was  usc-d 

for  each  solution.  Xo  convergence  problems  were  encountered. 

The  finite  difference  metliod  described  in  Section  IV  is  a  t  iru-  tU - 
pendent  technique.  The  solution  is  advanced  in  tint  by  I  i  nc  ri-mii  n  t  s  . 

Tile  numerical  solution  was  considered  to  have  reai-bed  a  stead-.-  state  wiai. 
chan.ges  in  the  computed  values  of  the  lift  and  drag  coe  t  f  i  t- i  eii  t  s  were  K  . 
then  1'  o\'cT  one  clia  rar  I  er  i  .s  t  i  c  or  nond  ime-n.s  i  ona  1  time  periis!.  I  icb  -  li.ad. 
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state  solution  required  three  to  six  characteristic  time  period.',  and  fro- 
two  to  four  hours  of  CPl'  time  on  a  CDC  CYBER  750  computer. 

The  numerical  solution  of  the  mean  flow  field  near  t'ne  airfoil  si  c- 
tion  at  each  angle  of  attack  is  presented  by  using  .streamline  oiintours, 
velocity  field  vectors,  and  pressure  contours.  The  mean  .''tream.l  jt-'.o  ci.mi- 
tours  are  shown  in  Figures  4  through  7.  A  small  laminar  separation  butoK 
on  the  upper  surface  with  negligible  trailing  edge  separation  is  seen  in 
Figure  4  for  a  =  5°.  For  an  increased  angle  of  attack  equal  to  7.5',  tl.e 
separation  bubble  has  decreased  in  size  and  moved  forward  toward  the  leading 
edge  of  the  airfoil.  Figure  5  also  shows  the  turbulent  trailing  edge  .'eparaiii: 
region  which  has  also  moved  forward.  Tlie  trailing  edge  separation  region 
has  grown  significantly  and  has  progressed  to  the  quarter  chord  pc\=ition 
when  the  angle  of  attack  has  reached  11.5'^.  This  behavior  is  typically 
observed  during  trailing  edge  stall  (72) . 

The  velocity  field  vector  plots  presented  in  Figures  8  through  11 
also  illustrate  the  phenomena  discussed  above.  In  addition,  a  sm; 1 1  scnlt- 
clockwise  rotational  motion  is  observed  in  tVie  rear  portion  of  the  trail  ir.g 
edge  separation  region.  A  circulatory  flow  pattern  witliin  the  lam.i’iar  supcira- 
tion  bubble  is  seen  for  =  5'’  in  Figure  8.  Ibc  boundary  layer  mar  the 
leading  edge  for  each  angle  of  attack  contains  at  least  five  grid  points 
along  each  ^  contour  while  approximately  ten  points  resolve  t!;e  lo-w(.r 
surface  laminar  boundary  layer  near  the  trailing  edge.  Thi.^i  numVu  r  .1  I'oir.i: 
should  be  sufficient  to  resolve  the  primary  features.  'lia  pres.san  >  o?:- 
tours  neat  the  airfoil  in  I'igures  12  through  15  clearly  show  i  xioi^.d.  i 
regions  of  low  pressure  above  the  airfoil  and  iiigh  prc.ssnri.  inline  tiie  .lir- 
foil  a.s  tli^-  ancle  of  attack  increases.  The  large  jiressuri  variation:  in 
the  upper  surface  suction  peaks  are  also  ob.served.  lin  [.ri  ssnri  i  it  lii 
indicate  that  only  sir, all  piressure  variations  occar  in  tin-  wain  tor  i  ' 
an,;le  o\  uttmk.  T!u  snoot!,  and  cor.timious  numeric.il  S'lnLit.''  lor  !  !a 


pressure  field  is  in  sharp  contrast  to  the  large  oscillatory  results 
obtained  by  Hodge  (7)  for  laminar  flow. 

For  a  moderately  thick  airfoil  with  a  smooth  surface  in  a  flow 
with  a  low  freestream  turbulence  intensity,  the  laminar  boundary  layer 
near  the  leading  edge  can  separate  upon  encountering  a  strong  adet-rst' 

■■  pressure  gradient.  Subsequently,  shear  layer  instabilities  can  cau^L 
transition  to  turbulence.  The  increased  mixing  may  reattach  the  shear 
layer  and  thereby  form  a  separation  bubble.  The  formation  of  separa¬ 
tion  hubbies  and  their  relationship  with  the  stalling  characteristics 
of  airfoils  at  various  Reynolds  numbers  have  been  investigated  by 
Gault  (85),  Caster  (86),  Hoad,  et  al  (87),  and  Arena  and  "ueller  (8f). 

The  laminar  separation  bubble  has  been  observed  experimentally  by  hregor;.' 
and  O'Reilly  (89)  for  the  KACA  0012  airfoil  over  a  large  range  of 
Reynolds  numbers.  These  experimental  and  empirical  results  will  be  used 
to  compare  witli  the  numerical  solutions. 

In  this  investigation,  the  beginning  of  turbulent  transition 

the  upper  surface  of  the  airfoil  and  the  transition  length  b.avc  been  ba -ec! 

on  the  closure  of  the  laminar  separating,  shear  layer.  The  .surface  mean 

pressure  coefficients  for  o  =  5“  and  7,5''  presented  in  i'igures  1  and 

17  clearlv  show  the  laminar  bubble  constant  pressure  re.- ion,  downs  t  n  a: . 

of  the  suction  peak,  followed  by  a  rapid  recovery.  Tran.sition  in  ti.r 

numerical  solution  occurs  near  the  down.stream  side  of  thu  bub’  I.,  w'ltie 

the  steep  pressure  recovery  begins.  This  phenomenon  ha?  hewn  r  p.  rtid  b-,- 

Wallis  (90)  and  Arena  and  Mueller  (88).  the  chord  location'  ler 

start  of  transition  and  full  turhulcnco  >:.|  rt-lative  to  thi  <  p.;  to  t  i  on 

and  rea  1 1  achiaen  t  points  of  the  bubble,  x  and  x,,  respt  t  t  i  vo  1  ,  an  gic.n 

S  r’ 

in  Table  1.  For  the  smaller  angles  of  attack,  rea  1 1  achr.c!;  t  o<.  nrttd 
prior  to  attaining,  fully  turbulent  flow.  This  behavior  ha?  lu  en  in- 
p(  r  i  men  t  a  1  1  obsorved  bv  .Aro:ia  and  Muoller  (8bi  iii  t  hi  it  1  ov.'  iovn.  !d: 

f.  7 
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TABLE  I 

Computed  Characteristics 

The  Airfoil  Leading  Edge 

a  (Deg) 

5.0 

1 

O 

ro 

7.5 

-  0.47 

9.5 

- 

11.5 

_ 

of  Turbulence  Near 
(x  =  -0. 5) 


X 

t 

0.24 

-  0.15 

-  0  ,  DO 

0.42 

-  0.38 

-  0.  35 

0.48 

- 

-  0.40 

0.49 

_ 

-  0.48 

TABLE  II  Conputed  Laminar  Separation  Bubble  Characteristics 


(Deg) 

“■I 

S 

»12S 

Ke .  . 

-  ^ 

■  s 

"R 

5 

.124 

3.1 

680 

.002  5 

.27 

-.002 

7.5 

.091 

3.2 

540 

]68 

.0013 

.09 

-  .002 

number  flow  study. 

The  coinputed  values  of  several  parameters  at  separation  for  both 

0=5°  and  a  =  7.5°  cases  agree  w’ith  empirical  results  and  are  presented 

in  Table  11.  The  pressure  gradient  parameter  M  =  -  '  Re  du  /t!s 

in  nond imensional  variables,  where  ’  is  the  momentum  thickness  and  s 

is  arclength,  compares  favorably  with  the  laminar  separation  criteria 

of  Curie  and  Skan  (91)  given  by  ^  0.09.  Tlie  shape  parameter  = 

i*/  takes  on  values  of  3.1  and  3.2  at  laminar  separation  compared  withi 

an  average  value  of  3.5  reported  by  Curie  and  Skan.  The  Reynolds  number 

at  separation  based  on  the  displacement  thickness  (Re;j.g)  lias  a  value 

greater  than  500  and  predicts  a  short  bubble  using  the  Owen-Klanfer  (^2) 

criteria  confirmed  by  Gault  (85)  and  Crabtree  (93).  Th.e  bubble  lengtii  1.^^ 

is  of  the  order  10“  C*g  and  decreases  rapidly  with  increased  angle  of  attack 

as  Seen  in  Table  II.  This  relationship  also  indicates  a  short  bubble  (90, 

U 

93)  whereas  long  bubbles  have  lengths  of  the  order  10  (85,8b).  The 

assumption  tliat  laminar  separation  precedes  turbulent  transition  is  verified 
with  the  computed  Reynolds  number  at  separation  based  on  momentum,  thickiuss 
Re  .  Crabtree's  (99)  criteria  indicates  that  transition  has  occurred  if 
Re  >  780  whc'n  first  reaches  0.09.  Thus,  the  solutions  correctlv  piudicl 
that  laminar  separation  occurs  before  transition  to  turbulence  .  e  turbu¬ 
lent  reattachment  criterion  for  the  bubble  given  by  Roberts.  (‘'5)  is  = 

(  ! w  )  du  /ds  >  -  0.0059.  This  criterion  is  alsci  satisfied  a,-.  bown  in 

e  e  - 

Talile  TI. 

The  trends  discussed  aliove  for  >i,. ,  ' ,  and  both  hi'undar-.-  laar  i.,  "nild;- 

S  In 

numbers  continue  at  the  higher  angles  of  attack.  Hewc  vi  r ,  at  ai:,;,l,c  el 
attack  near  stall  the  hublle  length,  decreases  ti'  about  1  t  iiord  (  'Oi  -..hie' 
is  the  ordi  1  of  till-'  ’  contour  spacing  neai  tlie  lead.iii;.'  id.i  .  '!  hi  iu:;i  r  i  >  s,  ' 
solution:  at  >  =  9.5  and  i  =  11.9  consequent  1  v  d.o  •so!  r.-o.-l.i-  to,  hoi’  ,  .. 


seer,  in  Fii;ure.s  18  and  19.  The  solutions  become  very  sensitive  to  small 
changes  in  the  turbulent  transition  length  and  starting  location  .  li  r 
example,  at  a  =  9.5'^  with  an  initial  transition  location  of  =  -('.48, 
a  change  in  transition  length  from  0.022  to  0.025  produced  shedding 
of  hubbies,  large  pressure  oscillations,  and  a  significant  30'  inerra,-, 
in  average  drag  compared  with,  the  steady  state  solution.  The  expcf  r  i::.L  r.  t  a  1 
pressure  data  were  measured  with  diaphragm  transducers  capable  of  measuring 
frequencies  well  beyond  5000  Hz.  No  dominant  frequency  was  observed. 

Hoad,  et  al  (87)  recently  observed  a  small  scale  unsteadiness  in  a  slurL 
region  near  the  surface  of  a  N.'\CA  0012  airfoil  at  angle  of  attack.  Ibe 
laser  velocinieter  mean  velocity  histograms  for  a  certain  region  near  tin 
leading  edge  exhibited  dominant  and  small  secondary  mean  velocities.  Ihcd, 
et  al  suggest  that  this  behavior  may  indicate-  an  unsteadiness  within  tin 
laminar  separation  bubble.  The  remaining  flow  field,  however,  was  steady. 

The  numerical  solutions  presented  for  each  angle  of  attack  are  tin- 
result  of  a  parametric  study  where  the  upper  surface  tranxition  length  and 
location  were  varied  in  order  to  close  the  bubble  and  maintain  .steady  f  hcv.- . 
In  eacii  case,  a  further  increase  in  tlie  transition  length  or  movenicn.t  down- 
.strear.  of  tin-  start  of  transition  resulted  in  a  non-ph.ys  ical  osc  i  1  1 t  -ry 
motion  propagated  down.stream  with  an  accompanied  large  increa.si  ii;  dr,;c. 

The  good  agreement  between  th.e  calculated  and  empirical  sL-p.i  ra  t  i  on  h,.’  1  li 
characteristics  indicates  that  the  origin  of  Lurbnience  o;.  tin  iqg.  r 
face  is  satisfactorily  ru'delled. 

Computed  boundary  layer  vi  1  oi  itv  profiles  at  four  slat  inn-  on. 
upper  surface  are  compared  with  two  sets  of  hot  wir.-  am  .’■.c>::.i.  t  r-,  da!,-  i  in 

Figures  fC  throu  'li  23.  ''he  profiles  were  measured  an.d  iorqutid  iit  Ih.  -  a- < 
locations  on  const.ant  lines  loiaited  .approx  imat  i  1 '■  at  18  ,  'ii  -  ,  f  ,  .ind 

“2  rhord.  I’hi'  pr.'fileo  are  nond  i  men  s  i  on.; !  i  ced  !■■■  the  h-iMllg  oi-'  n;:li  ' 
houiidarn  la-.a-r  1 1.  i  okm  -o-  which  is  -iefim-d  a,-,  th.  n-orn.il  .ii'-tim.i  ir-:-  ib. 

7  <  ■ 


iirTftii'  <1  .  i  n  .j 
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surface  to  where  the  tangential  velocity  f  ir.  t  attiint  '-‘'-j  e:  ti..,  ’  . 

edge  velocitv.  Significant  Kf.S  fluctuations  for  eae;.  a:;gii  <  :  ati.;.  .  w. 
detected  experimentally  at  the  inner  data  point  ]i)cetJ<e,  :<  r  tati.-c  is 
through  four  whic’”.  indicate  the  presenci.-  of  a  tur'-sIeM  uc-dary  hv'  i  . 

No  fluctuations  v.-ere  ohse-rved  in  thn'  lower  surlaci  iw  uc,:  a  r ;  lie,,  r  .;t 

t'nt  9'r  ctiord  location.  This  re.sult  provides  e.xperi::...  ;.ta  i  tcici  ;  i.  :  ; 

the  assumption  t’nat  the  lower  surface  boundary  la;.'er  ris.aixs  in 

the  presence  cf  the  favorable  pressure  gradient  at  the  -  i  ;--  Ids  nu:d  i  i 

under  investigation.  The  calculated  boundar;.'  layer  bt-comt-s  thirinr  dew:.- 
stream  compared  with  the  experimental  profile.s.  li.is  te.-clt  iwcIl:  it 
caused  by  grid  boundary  layer  resolutitu.  i-rrors  pro:a,'a;,,d  c.ow:;.- 1  re.g  .r 
by  deficiencies  in  the  tubulence  modul.  The  iiOt  wirt-  d.ita  iia'.a  a:;  t  t  ^ - 
mated  experimental  error  of  5  . 

The  boundary  layers  on  the  upper  and  lower  surfaee.s  merge  to  form  tin 
wake  downstream  of  the  airfoil.  Wake  rt.ea:;  veloi  ity  profiles  lU  chcimi  lo¬ 
cations  of  0.58,  0.79,  and  1.54  for  c'ach  ancK  ef  attarV  art.  shown  ir. 
Figures  24  through  27.  Tiie  trailing  edgt^  i.s  locati-.!  at  x  -  i  .S.  pv  - 

file.s  are  measured  along  constant  '  I  ir.e.s  with  the  t'rigir.  on  tht.  extindt  d 
chordline.  Th.e  tSticker  numerical  upper  .surface  houndary  larcr  r  iipa  r e;- 
downstream.  Tiie  computed  large  variation.?  of  velocity  near  lin.^  h.'Wt  i  tc.  t 
of  the-  wake,  whicli  come  from  the  lower  .surface  boundary  lavt  r,  art  in.  no 
av,reer;ent  with  the  hot  wire  anenometry  measurements.  the  mi-as  ■ .  r  er  ,  . ,  t  -  vaia 
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i  Kr'Si'IJr,  stri'^s  u  v  ,  noncira-n^  iua.i  1  i  i.'saalij 

aad  (.asitour  i)lal-  ar^'  shovn  in  1‘isiiros  Ad  Ihrourk  H  -u-  SL-o-ia  iri. 
jKiltnr:;  i;^  qua  1  i  ta  t  1  v  u  1  >  similar  to  t  tit'  oxpcrir.icnta  1  data  cdataincA  ay 
Colo.s  arai  v.'adaock  For  a  I.A'vCA  k4i;;  airfoil  waicii  art-  siu'v,"-.  is 

i'iyarf  'AI  .  Sor.t-  exj'or  irifn  t  al  t:t-asurer:fnts  v.-ero  cade  durir.y  tho  ;v.'i 
vire  anor:OA.oCry  study  (84)  and  are  presented  in  Appendix  '!:.;s  d/Aa 

also  illustrates  the  rapid  variation  of  negative  and  positive  stressi- 
in  the  upper  and  lower  portions  of  the  near  wake,  resi)eft  ive  ly .  .'A-.-:.:- 

tudes  of  tlie  order  0.01  have  been  observed  (84,  96)  in  the  near  wake  and 
Were  also  obtained  numerically  where  the  contour  values  vary  fro:v.  -('.(A  li 
+  ■'.01  in  Figures  28  tliroug'n  11. 

The  comqiutfc'd  surface  niean.  pressure  distributions  for  all  four  cinyli, 
of  attack,  presented  in  Figures  I''  throug'i  !  avorahly  compart  wit': 

tVie  ex'perinen  tal  data  (84)  which  have  an  esticated  t:rror  of  1  to  f  .  A", 

invisoid  solution  with  imposed  Kutta  condition  wa'-'  computed  usiiyy  t:.e 
approaoh  outlined  in  ’.ppendix  G  and  is  also  shown.  The  com.  ..t  ed.  Aa'.'irr- 
.‘fto'kes  lar.inar  suction  peak  is  well  defined  and  tncreasf,s  rapid:;,  wit:, 
ancle  of  al'ack.  The  lower  surface  pressure  pe.ik  !iot!.  r.o.'e.-  do'wr  tri.i:' 
and  liecoc.t  ■  brtiader  as  tine  angle  c'f  attack  incre.ises.  I  :  is  iu  :.  :.  :-  r  A 
also  seLv:  in  tiie  exjuT  Imen  ta  1  data.  A  separ.ition  hu''lt  i  o’  .-i  rv- 

the  ex['e  r  i  nv  II  ‘  .1 1  d.ata  where  the  curvature  rapiiilv  etian,;.!  .  .:  '  Iv.c: 

occurs  furtiier  downstr<'ar  when  eom;[:,ired  wit'n  the  Awcier-AA  s.  ic  ju', 

Till-  dif'Crencc  is  prid'ably  caused,  lyv  f  re-e  - 1  it  a:'  I 'i  ri-u  1 1  i:  ■  t  w!.;-  i.  '  vl' 
i('::  an'.;  ti.i-  .t uAi-.t-qiit 'U t  fcu"‘:atio:'  c-f  the-  l’u:''l-*'-  .  •'  t  *.  '  i  ,  • 
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data  vv’erc  nieasured  in  a  wind  tunnel  vit't:  a  freesLreat.  t  ur  l-u  i  enc  e  inte:..-  il 
o:  approximately  0,5  compared  with  t’le  unperturbed  freestrear.  of  t'ae  r.i;'  - 
erical  solution.  Some  disagreement  occurs  in  the  trailing  edge  region 
whore  the  computed  result  ha.s  a  more  pronounced  flat  trailing  edge  stall 
character  i.st  U'  . 

The  computed  inviscid  surface  pressure  d  is  t  r  ibul  io:^  cons  i  s  t  exi  L  1  ■■ 
predicts  both  a  larger  suction  peak  on  the  upper  surface  .and  a  largor  pr.,  s- 
sure  peak  on  the  lower  surface.  The  discrepancie.s  increase  vcith.  angle  of 
attack.  The  inviscid  solution.s  predict  complete  recovery  at  the  trail  in. 
edge  wit'n  a  stagnation  point  and  =  1.0.  The  large  pre.s.surg-  grarl  ion.  L 
in  the  inviscid  and  viscous  solutions  near  the  airfoil  r.o.se  are  rt.solvcd 
wh.ich.  indicates  an  adequate  grid  point  distribution  in  the  regio:-.. 

.-\n  exam.in.at  ion  of  the  experimental  and  inviscid  .surface  pros.-uro 
d  istr  ibut  ions  for  zero  angle  of  attack  given  in  Figure-  33  provides 
additional  infonrg.tion.  The  inviscid  pressure  distribution.s ,  obtained 
withi  the  finite  alfference  techni'':ui*  in  .Appendix  C,  1‘or  tlie  upper 

end  loxL-r  .surfaces  are  indistinguishable.  This  physically  corre.t 
syr.m.L.t  r  ic  result  cam.v'  directly  iron  the  numerical  solution,  and  net 

.specified  a.s  a  boundary  condition.  Fu  r  thermo  re ,  the  num.erical  in.\'is.'id 
result  is  virtually  identical  with,  an  inviscid  solution  ('’7'>  u.sing  t 'n.e 
method  of  Theodorsen  (12).  The  experimental  data  for  hot:,  surfaci:- 
agree  well  with  the  inviscid  result.s  wh.ich  indicates  the  sm..;ll  effe.l 
tlnit  V  i  ■  ,'s  i  t  y  ti.is  on  the  pre,-.sure  distribution  for  ,i  s  t  ri  am  1  i  t'.i  - 

metric  airf'.'il  ut  zero  angli,  of  attack.  This  agrcenteTtt  al.  j  rovid.. 
evidertre  that  ttie  expo  r  imeit  t  .i  1  airr.-'il  section  dots  ap;  ■■oxir  a  ;  ;  ygo' 

Onlf  mn',  i.^iu.itiou.  V'iscou.'  iffert‘-  c'n  tVii-  pressr.rt  d  i  t  r  i !  u '  h-n  d.n  .  ..or 
(.IS  W.1-  p  r  o' '  1  ni  I  ■  I  ■■  nentioie.!'  .is  the  artzli-  of  attaoK  is  intio.a-i  !  i  -.'on  t.o 
‘t'.o’l  .  .  omjatis.'n  hetwee::  the  ospe  r  i  r.e!- 1  a  1  dat.i  ;...r  th.  tw  ‘o;,-.  ,  , 


f 
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rt'X'eals  a  si  ig!u  suction  pt-ak  on  1)10  surface.  '!  i.t-  suction,  jsak 

sugi'osts  that  tlie  experinontal  zero  angle  of  attack  may  actual  1;.'  l-e 
slightly  negative.  A  snal  deviation  in  anvle  is  plausible  since  t!;e 
expt'r ir.en t a  1  zero  was  de t ert.iined  by  manually  ad’ucti!c.;  the  airfoil 
attitude  until  the  pressure  distributions  c^n  loth  surface.-^  apnie.r-. 
identicc;!  . 

The  lift  cc'ef  f  ic  ients  obtained  f  roiv.  thie  present  turbul-.-nt  he.v  i  i  !' - 
StokciS  solution,  int'iscid  results,  and  two  sets  of  exper imen tc;  1  data  arc 
cot:pared  in  Figurc^  3A.  The  experimental  results  sh;"wn  have  In-cn  c -c  rre.,  t  ■.  .. 
for  angl^j  cf  attack.  The  effective  experimental  zero  ancle  c  f  attecl.  f'r 
a  symc'-et  r  ic  airfoil  section  is  defined  to  be  t'c.e  ancle  where  the  1  i t 
coefficient  is  .zero.  Fi-r  th.e  Seiler  l.aborator;.’  results  the  c:f.,c- 

tive  zero  was  found  to  re  0.5'.  This  snail  ccrrectiox  is  in  accerd  with 
the  analc'sis  of  the  nominal  zero  angle  of  attack  pre.ssure  d  i  s  t  r  i  bu  t  i  on 


results  already  pre-'sent  ed  . 


xperinc-r. tal  data  have  t'-.c-ref c-i'e  V 


translated  0.5"  in  tl.e  necative  ■  direction  in  I'icure  '-s.  The  h.’.CA 


results  (.s'-A  were  sir.ilarlv  nubllshc-'d  in  corrected  fore 


'tbiOr  rc~-en 


twci-d  ir.ens  1  ona  1  wind  tunnel  corrections  including  solid  biockinc,  wake' 
Mocking,  and  streamline  curvature  effects  Were  investigated.  h:-.;- :  r  i  c,M 
low-speed  wind  tunnel  Ccorrection  fartc'rs  (‘■'7!  using  the  Sc  ilc-r  1  h'  r,:'.  r 
d'x3'  tunnc'l  geometry  were  applied  to  the  expor ic.en t a  1  data  at  11.' 
account  for  blockage  and  wall  i  n  t  er  f  >  rer.ce  .  ’Ihe  corri- t  c-d  Vj'uc  •  : 


the  lift  cot-f  f  ic  lent  and  an.cl. 


lift  c' Oi' f  f  i  r  i  en  t  and  +  A.’  ,  r- 


Ficure  have  been  so  correct- 


L 


’f  r  ir.'.t 


NA'.',.\  Ai*  :  (  ^  i".  l;i  i  a'l 
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The  NACA  lift  curve  then  flattens  more  rapidly  near  stall.  These  minor 
differences  may  be  explained  by  a  larger  freestream  turbulence  level  in 
the  NACA  results  which  tends  to  delay  the  turbulent  trailing  edge  separa¬ 
tion.  Also,  the  major  variation  of  the  lift  curve  at  small  angles  of 
attack  occurs  for  Reynolds  numbers  less  than  10^  (56,  89).  Freestream 
turbulence  can  cause  more  rapid  transition  and  result  in  an  apparent 
freestream  Reynolds  number  greater  than  the  nominal  turbulence-free  value. 

The  computed  turbulent  Navier-Stokes  lift  coefficients  for  the 
four  angle  of  attack  cases  are  in  excellent  agreement  with  the  experi¬ 
mental  data  as  shown  in  Figure  34.  The  force  components  were  computed 
using  trapezoidal  integration  to  sum  the  total  surface  stresses  obtained 
from  the  flow  field  solution.  The  numerical  results  exhibit  a  gradual 
decrease  in  slope  with  increased  angle  of  attack  similar  to  the  Seiler 
Laboratory  data.  This  behavior  is  consistent  with  the  data  because  the 
numerical  solution  has  a  quiescent  incoming  freestream  with  a  correspond¬ 
ing  smaller  apparent  Reynolds  number.  The  significant  effect  that  the 
viscous  separation  in  the  trailing  edge  region  has  on  the  lift  coefficient 
for  angles  of  attack  greater  than  8°  is  observed.  The  computed  values 
are  within  5%  of  the  experimental  results.  The  lift  coefficient  was  also 
computed  for  each  angle  of  attack  using  the  contour  momentum  integral 
method  described  in  Appendix  D.  The  calculated  values  differ  by  less 
than  1%  for  ri  contour  paths  within  one-half  chord  of  the  aifoil  surface. 
This  result  demonstrates  the  near  flow  field  consistency  in  resolving  the 
lift  force.  The  time  dependent  terms  in  Equations  83  and  84  always  con¬ 
tributed  less  than  0.5%  of  the  total  computed  lift  which  indicates  again 
a  converged  steady  solution. 

Two  inviscid  flow  predictions  for  the  lift  coefficient  are  shown  in 
Figure  34  for  comparison.  The  numerical  inviscid  finite  difference  results 
were  computed  with  the  second  order  accurate  trapezoidal  integration 
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technique  used  in  the  viscous  calculations.  The  linear  airfoil  theory 
2tt  slope  is  included.  The  inviscid  flow  results  over-estimate  lift  by 
25%  at  an  angle  of  attack  of  5°.  This  discrepancy  increases  with  angle 
of  attack  because  the  viscous  effects  which  induce  stall  are  not  pre¬ 
sent  in  an  Inviscid  calculation. 

The  computed  drag  coefficients  for  the  four  cases  are  compared  with 
available  experimental  data  (56)  in  Figure  35.  The  force  components  were 
computed  using  trapezoidal  integration  to  sum  the  total  surface  stresses 
obtained  from  the  flow  field  solution.  The  drag  component  is  in  the 
direction  of  the  incoming  freestrear;  flow.  The  rapid  increase  in  the 
drag  coefficient  which  accompanies  the  smaller  increase  in  lift  at  higher 
angles  of  attack  near  stall  is  observed.  Jacobs  and  Sherman  (56)  acknow¬ 
ledged  that  the  experimental  drag  data  in  this  Reynolds  number  region 
have  error  greater  than  the  estimated  +  0.001  for  the  larger  Reynolds 
number  results.  The  agreement  between  the  present  numerical  solution  and 
the  experimental  data  is  within  ten  drag  counts  in  the  region  of  the 
maximum  lift  to  drag  ratio. 

The  presented  numerical  solutions  for  two-dimensional  incompressible 
turbulent  viscous  flow  over  airfoils  were  obtained  with  several  parametric 
studies.  The  effects  associated  with  varying  the  turbulence  transition 
values  and  the  SOR  iteration  parameters  including  time  step  size  have 
already  been  discussed.  The  influence  of  other  turbulence  model  para¬ 
meters,  far  field  boundary  conditions,  and  grid  fineness  on  the  numerical 
solutions  have  also  been  investigated. 

The  turbulence  model  contains  several  parameters  which  can  vary.  The 
effects  that  changes  in  these  parameters  have  on  the  flow  field  near  the  air¬ 
foil  surface  and  on  the  values  of  lift  and  drag  were  examined.  The  value 
of  the  nominal  inner  eddy  viscosity  parameter  k^^  was  increased  by  20°'  at 
a  5°  angle  of  attack.  The  lift  coefficient  subsequently  increased  by  2% 
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while  the  drag  coefficient  increased  by  57  or  six  drag  counts.  These 
integrated  effects  are  caused  by  the  observed  small  increase  of  velocity 
in  the  inner  boundary  layer  on  the  upper  surface.  The  outer  eddy  vis¬ 
cosity  parameter  was  next  increased  by  50%.  The  result  was  an  increase 
in  turbulent  mixing  with  a  corresponding  thicker  boundary  layer.  The 
velocities  near  the  edge  of  the  upper  surface  boundary  layer  decreased 
which  caused  an  87  loss  of  lift  with  no  apparent  effect  on  the  calculated 
drag. 

Values  for  the  constants  in  the  inner  eddy  viscosity  parameter  k^ 
relaxation  factor  given  by  Equation  65  were  obtained  from  a  parametric 
study  at  an  angle  of  attack  of  11.5°.  This  angle  of  attack  was  chosen 
because  the  effects  of  the  relaxation  are  significant  only  for  angles  of 
attack  near  stall.  A  solution  was  initially  computed  without  the  relaxa¬ 
tion  factor  f.  The  laminar  separation  bubble  remained  closed  with 
trailing  edge  separation  beginning  at  mid-chord.  The  calculated  values 
for  the  lift  and  drag  coefficients  were  0.98  and  0.064,  respectively.  An 
examination  of  the  lift  curve  in  Figure  34  reveals  that  this  solution  ex¬ 
hibits  the  character  of  leading  edge  stall  where  an  approximate  linear 
behavior  is  sustained  until  separation  abruptly  occurs.  Leading  edge  stall 
does  occur  for  the  NACA  0012  airfoil  for  Reynolds  numbers  greater  than 
about  500,000  (56,  89).  Thus  the  mechanism  for  leading  edge  stall  seems 
present  within  the  modified  turbulence  model.  A  numerical  solution  was 
next  obtained  using  the  relaxation  factor  with  a  relaxation  distance  s^  = 
0.25  and  a  delay  distance  s^  =  0.  in  nondimensional  distances.  The  laminar 
bubble  remain  closed,  but  the  trailing  edge  separation  region  was  signifi¬ 
cantly  larger  and  moved  to  within  the  20%  chord  location  (x  =  -  0.3).  The 
lift  coefficient  decreased  substantially  to  a  value  of  0.8  while  the  drag 
coefficient  increased  moderately  to  a  value  of  0.074.  The  final  solution 
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which  has  been  previously  discussed  was  calculated  using  distances  = 

0.5  and  s^  =  0.1.  The  computed  lift  and  drag  coefficients  became  0.84 
and  0.068,  respectively,  and  are  in  agreement  with  the  experimental  data. 
The  sensitivity  of  the  numerical  solution  to  these  parameters  has  thus 
been  obtained.  The  last  set  of  distance  parameters  was  used  in  calcu¬ 
lating  the  solutions  for  the  other  angles  of  attack.  The  onset  of 
trailing  edge  separation  for  the  5°  and  7.5°  solutions,  seen  in  Figures 
4  and  5,  indicates  that  the  relaxation  factor  was  not  activate.d.  Thus 
the  relaxation  model  for  parameter  affects  the  flow  field  only  near 
stall  where  the  pressure  gradient  approaches  zero  over  a  large  portion 
of  the  upper  surface. 

The  relaxation  distance  in  the  near  wake  turbulence  model  given 
by  Equation  74  was  also  investigated.  Numerical  solutions  were  obtained 
at  5°  angle  of  attack  for  relaxation  distances  equal  to  5  and  100  times 
the  lower  surface  boundary  layer  thickness  near  the  trailing  edge.  No 
changes  in  either  the  surface  pressure  distribution  or  the  integrated 
force  coefficients  were  observed.  In  two  recent  near  wake  calculations 
with  similar  relaxation  models,  Waskiewicz  (98)  used  a  value  of  30  trail¬ 
ing  edge  boundary  layer  thicknesses  and  Hasen  (99)  used  10  boundary  layer 
thicknesses.  The  relaxation  distance  of  5  boundary  layer  thicknesses  was 
retained  since  the  agreement  of  the  Reynolds  stress  field  discussed  pre¬ 
viously  was  improved  for  this  choice. 

The  effect  on  the  solution  of  the  placement  and  type  of  far  field 
boundary  conditions  should  be  considered  in  any  computational  work.  In 
this  investigation,  two  types  of  boundary  conditions  at  different  loca¬ 
tions  were  examined.  The  four  angle  of  attack  solutions  that  have  been 
discussed  were  computed  using  the  freestream  far  field  boundary  conditions 
described  in  Sections  III.B  and  IV. C  at  a  circular  outer  boundary  of  radius 
10  chords.  The  effect  of  the  far  field  boundary  placement  on  an  Inviscld 


solution  for  Joukowski  airfoils  was  previously  accomplished  (82) . 

This  analysis  indicated  that  a  radius  of  10  chords  was  sufficient  as 
discussed  in  Section  IV. C.  The  effect  of  the  boundary  placement  on 
the  present  viscous  solution  was  examined  by  using  an  alternate  outer 
boundary.  The  J  =  40  n  contour  with  a  semi-major  x  axis  of  4,77  and 
semi-minor  y  axis  of  4.64  was  chosen  as  the  new  outer  boundary  because 
it  approximates  a  circle  with  half  the  previous  radius.  The  numerical 
solution  was  then  obtained  for  ot  =  5°  and  compared  with  the  solution 
using  the  10  chord  radius  outer  boundary.  The  average  values  for  the 
lift  and  drag  coefficients  were  identical.  The  solution  with  the  closer 
far  field  boundary  had  a  small  oscillatory  behavior  with  variations  in 
the  lift  and  drag  coefficients  of  1%  and  2%,  respectively.  Further  ex¬ 
amination  of  the  near  surface  flow  field  revealed  that  the  laminar 
separation  bubble  had  a  small  scale  unsteadiness  which  locally  affected 
the  pressure  field. 

The  far  field  potential  flow  boundary  condition  model  developed  in 
Appendix  F  was  next  applied  with  the  outer  boundary  locations  of  10  chords 
and  about  5  chords  in  turn.  This  model  approximates  the  small  pertur¬ 
bations  from  the  freestream  conditions  for  the  velocity  and  pressure 
fields  at  a  large  but  finite  distance  caused  by  the  presence  of  a  body  in 
the  flow  field.  In  this  approach,  the  upstream  boundary  conditions  for 
both  velocity  and  pressure  became  the  calculated  values  from  the  inviscid 
potential  flow  model.  The  downstream  boundary  condition  for  pressure 
also  became  the  corrected  value  rather  than  the  freesteam  value.  The  no 
change  downstream  boundary  condition  for  the  velocity  components  was  re¬ 
tained.  The  solution  was  calculated  with  the  revised  far  field  boundary 
conditions  at  the  10  chord  radius  circular  outer  boundary.  The  boundary 
conditions  were  initially  updated  every  200  time  steps  (0.2T)  by  using 
the  latest  calculated  value  of  the  circulation.  After  three  characteristic 
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time  periods  T,  the  mean  values  of  the  lift  and  drag  coefficients 
were  '/irtually  identical  with  the  previous  freestream  boundary 
condition  results.  However,  while  the  lift  coefficient  had  variations 
of  less  than  1%  from  the  mean  value,  the  calculated  drag  coefficient  was 
periodic  with  a  period  of  0.4T  and  variations  of  +  10^:.  Since  the  period 
of  the  oscillations  was  twice  the  period  of  the  boundary  condition  update 
and  the  circulation  was  essentially  a  constant  at  this  time,  the  outer 
boundary  conditions  (except  the  no  change  downstream  condition)  were  kept 
at  the  current  values  and  the  solution  was  advanced  1,5  characteristic 
time  periods.  At  this  time,  the  lift  coefficient  had  a  steady  value  of 
0.425  and  the  drag  coefficient  became  0.011  +  2%  (two  drag  counts). 

These  computed  results  are  within  12-  of  the  corresponding  coefficients 
calculated  from  the  solution  using  the  simple  freestream  conditions. 

Thus  the  more  accurate  boundary  condition  which  includes  the  effect  of 
a  finite  distance  from  the  body  induces  very  small  changes  in  tlie  computed 
flow  field  near  the  airfoil. 

The  modified  far  field  potential  flow  boundary  conditions  were  also 
used  at  the  near-circular  five  chord  radius  outer  boundary  defined  above 
to  obtain  a  solution  again  for  a  =  5° .  After  two  characteristic  times, 
the  outer  modified  freestream  boundary  conditions  were  held  constant  since 
the  lift  coefficient  had  small  variations  less  than  1%.  The  solution  ex¬ 
hibited  a  small  oscillatory  behavior  similar  to  the  freestream  boundary 
condition  result  even  after  four  characteristic  time  periods.  The  lift 
coefficient  had  a  mean  value  within  1%  of  the  freestream  boundary  condition 
value  with  variations  below  1%.  The  drag  coefficient  had  the  same  average 
value  as  the  previous  result,  but  the  oscillations  were  larger  with  a  de¬ 
viation  from  the  mean  of  approximately  +  10".'  (10  drag  counts).  The  un¬ 
steadiness  was  again  observed  to  emanate  from  the  laminar  separation  bubble. 


Small  pressure  oscillations  propagated  downstream  on  the  upper  surface 
of  the  airfoil  and  were  eventually  damped  out  near  the  trailing  edge.  In 
eac 1  case  no  changes  were  made  in  any  of  the  turbulence  or  other  parameters 
during  the  computations.  Small  changes  in  the  turbulence  transition  near 
the  bubble  would  probably  prevent  the  unsteadiness. 

The  far  field  boundary  condition  study  has  indicated  that  the  use 
of  the  freestream  boundary  conditions  at  the  large  but  finite  distance 
from  the  airfoil  surface  is  sufficient  for  the  computation  of  the  near 
surface  flow  field  and  force  coefficients.  Variations  in  the  outer  bound¬ 
ary  placement  and  the  use  of  a  more  accurate  far  field  boundary  condition 
had  negligible  effects  on  the  present  numerical  solution. 

The  numerical  implementation  of  the  downstream  no  change  boundary 
condition  was  also  investigated.  Second  order  accurate  central  spatial 
and  upwind  differences  were  used  in  separate  solutions  for  an  angle  of 
attack  of  7.5°.  No  difference  (+0.01%)  was  observed  in  the  surface  pres¬ 
sure  distributions  or  computed  force  coefficients  between  the  two  bound¬ 
ary  conditions.  The  computed  velocities  at  locations  across  the  wake  on 
the  outer  boundary  for  the  two  cases  differed  by  less  than  1%.  The  upwind 
difference  formulation  was  chosen  because  of  the  associated  convective 
properties . 

The  sensitivity  of  the  numerical  solution  to  the  fineness  of  the 
grid  was  examined.  The  7.5°  angle  of  attack  case  was  chosen  because 
both  the  laminar  separation  bubble  and  a  region  of  separated  flow  near 
the  trailing  edge  are  present.  A  coarse  grid  was  obtained  from  the  79x44 
grid  by  deleting  the  odd  numbered  n  contours  except  for  the  body  contour. 
The  coarse  79x23  grid  was  a  subset  of  the  full  coordinate  system  and  was 
used  to  study  the  effect  of  boundary  layer  resolution  on  the  numerical 
solution.  The  turbulence  parameters  and  SOR  acceleration  parameters  were 
unchanged  from  the  previous  solution.  The  coarse  grid  computation  w.-is 
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carried  out  over  four  chacteristic  time  periods.  The  calculated  values 
for  the  lift  and  drag  coefficients  approached  a  "steady"  periodic  state 
after  two  chaiacterisc  times  with  a  period  of  0.6T.  The  mean  value  of  the 
lift  coefficient,  averaged  over  the  last  two  characteristic  times,  became 
0.606  +  2/'  compared  with  the  full  grid  solution  value  of  0.606.  The 
computed  mean  drag  coefficient  was  0.0260  +  15%  compared  with  the  79x44 
grid  solution  result  of  0.0240.  Values  for  the  wake  velocities  in  the 
freestream  direction  along  the  outer  boundary  for  each  solution  compared 
within  1%.  The  unsteadiness  was  found  to  originate  from  the  oscillating 
laminar  separating  shear  layer  which  forms  the  front  portion  of  the  bubble. 
The  edge  of  the  boundary  layer  varied  between  the  J=3  and  J=4  r,  contours, 
as  defined  by  the  79x23  grid  system,  at  the  1%  chord  location.  This  motion 
also  caused  the  bubble  to  vary  in  length  and  produced  pressure  oscillations 
along  the  upper  surface.  This  unsteadiness  is  probably  attributed  to  the 
reduced  resolution  of  the  bubble  in  the  ri  direction.  Similar  sensitivity 
has  already  been  discussed  concerning  the  resolution  of  the  extremely 
short  bubble  in  the  f  direction  at  larger  angles  of  attack. 

The  computed  flow  field  near  the  airfoil  and  the  calculation  of  the 
force  coefficients  using  the  coarse  grid  are  in  agreement  with  the  results 
for  the  complete  79x44  grid.  A  quadratic  Richardson's  extrapolation  on  the 
computed  drag  coefficients  indicates  an  error  of  +0.0007.  This  agreement 
indicates  that  the  solution  obtained  with  the  full  79x44  coordinate  system 
is  adequate  to  obtain  +10  drag  counts  and  therefore  sufficiently  approximates 
Che  limiting  numerical  solution  obtained  from  an  exteraely  fine  grid.  Tlie 
coarse  grid  computation  further  indicates  that  the  numerical  multi-grid 
approach  (100)  may  be  Implemented  for  complex  flows  using  general  grid  trans¬ 
formations.  The  use  of  only  one-half  the  grid  points,  for  instance, 
during  most  of  the  time  marching  procedure  would  significantly  reduce 
the  computer  time  required  for  an  equivalent  numerical  solution. 


SECTION’  VI 


CONCLUSIO  .’S  AND  RECOMMENDATIONS 

Numerical  solutions  have  been  obtained  for  two-dimensional  incom¬ 
pressible  turbulent  flow  over  airfoils  near  stall.  The  time  dependent 
Reynolds  averaged  Incompressible  Navier-Stokes  equations  in  the  prim.it ive 
variables  of  velocity  and  pressure  together  with  a  Poisson  pressure 
equation  are  numerically  solved.  An  algebraic  eddy  viscosity  approach 
is  modified  for  separated  adverse  pressure  gradient  flows  and  used  to 
model  turbulent  closure  of  the  laminar  separation  bubble  and  the  subsequent 
turbulent  boundary  layer.  A  deficiency  in  the  standard  model  is  detected 
and  corrected  by  using  a  "limiting"  technique.  A  body-fitted  coordinate 
system  is  numerically  transformed  to  a  rectangular  grid  in  the  computational 
plane.  The  set  of  transformed  partial  differential  equations  is  solved 
with  an  implicit  finite  difference  method.  Successive-over-relaxation 
iteration  is  us<=d  to  solve  the  system  of  linearized  difference  equations 
at  each  time  step. 

Numerical  solutions  are  presented  for  a  N.ACA  0012  airfoil  near  stall 
at  a  chord  Reynolds  number  of  170,000.  A  sliort  laminar  separation  bubble 
located  near  the  upper  surface  suction  pressure  peak  is  obtained.  Computed 
laminar  separation  bubble  characteristics  including  the  criteria  for  sep¬ 
aration,  bubble  type,  and  turbulent  transition  agree  with  empirical  re.sults. 
Surface  mean  pressure  distributions  are  presented  and  found  to  comiiare 
favorably  with  experimental  data.  The  separation  bubble  is  observed  for 
angles  of  attack  of  5°  and  7.5°.  For  larger  angles  of  attack,  the  small 
bubble  essentially  disappeared  within  the  numerical  resolution  of  the 
streamwise  grid  spacing.  The  steep  leading  edge  suction  pressure  peak  is 
well  defined  for  each  angle  of  attack.  Velocity  profiles  at  four  stations 
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along  the  upper  airfoil  surface  are  compared  with  experimental  results. 

The  experimental  data  were  obtained  at  similar  grid  point  locations 
so  that  interpolation  was  not  required.  The  calculated  lift  and  drag 
coefficients  are  in  excellent  agreement  with  the  experimental  data. 

The  lift  coefficients  are  within  5"'  of  the  experimental  values  near  stall, 
and  the  computed  drag  coefficients  are  within  10  drag  counts  in  tlie  region 
of  maximum  lift  to  drag  ratio.  The  effects  of  viscous  separation  on  the 
lift  coefficient  curve  which  produces  a  maximum  value  are  seen  and  con¬ 
trasted  with  a  numerically  obtained  inviscid  potential  solution  with  Kutta 
condition.  The  observed  phenomena  of  trailing  edge  stall  is  predicted 
where  the  rear  separation  point  moves  forward  with  increasing  angle  of 
attack . 

The  sensitivity  of  the  numerical  solution  has  been  examined  in 
several  areas.  A  far  field  potential  flow  boundary  condition  which 
modifies  the  freestream  conditions  for  use  at  a  large  but  finite 
distance  from  the  airfoil  was  considered  for  two  outer  boundary  placements. 
The  study  showed  that  the  use  of  the  infinite  freestream  conditions  at  an 
outer  circular  boundary  of  radius  10  chords  produces  negligible  influence 
on  the  near  flow  field  and  calculated  values  for  the  force  coefficients. 

A  coarse  79x23  grid,  compared  with  the  79x94  grid,  was  used  to  evaluate 
truncation  error.  An  analysis  of  convergence  for  successive-over¬ 
relaxation  iteration  predicts  an  upper  limit  on  the  time  increnn-nt  for 
the  implicit  finite  difference  method.  Numerical  experiments  confirmed 
this  upper  bovind.  Several  parameters  within  the  turbulence  model  were 
systematically  varied.  The  only  significant  sensitivity  occurred 
near  the  downstream  side  of  the  laminar  separation  bubble.  Small  changes 
in  turbulent  transition  length  and  location  were  found  to  s Igni f ican t 1 > 
change  the  flow  field.  This  sensitivity  may  he  related  to  incipient 
turbulent  separation  which  results  in  failure  to  close  the  separation 
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bubble.  This  leading  edge  stall  phenomena  is  observed  at  higher 


Reynolds  numbers  for  the  NACA  0012  airfoil. 

The  agreement  between  the  numerical  solu’ ions  and  the  experimental 
data  and  empirical  results  indicates  that  the  eddy  viscosity  approach 
modified  for  separated  adverse  pressure  gradient  flows  adequately 
models  the  turbulence  on  the  upper  airfoil  surface.  The  near  stall 
airfoil  aerodynamic  characteristics  are  consequently  accurately  pre¬ 
dicted  by  the  numerical  method. 

The  results  of  this  investigation  have  suggested  areas  for  further 
research.  A  more  exhaustive  study  of  the  laminar  separation  bubble 
should  be  accomplished  to  better  understand  the  phenomena  of  leading  edge 
stall.  The  computation  of  a  geometry  with  detailed  experimental  results 
for  this  purpose  is  desirable.  Further  research  in  adapting  the  multi-  1 

grid  technique  would  significantly  increase  the  computational  efficiency 
and  appears  feasible  as  a  result  of  the  coarse  grid  calculation.  Turbulence 
modelling  is  an  area  of  substantial  interest  and  should  be  pursued  in 
conjunction  with  experimental  investigations.  The  improved  accuracy  of 
current  computational  aerodynamics  also  indicates  that  the  errors  in  ex¬ 
perimental  measurements  should  be  analyzed  in  a  more  rigorous  manner  so 
that  comparisons  can  be  better  evaluated.  The  method  should  also  he 
employed  in  time  dependent  problems  to  further  exploit  the  time  marching 
technique . 
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FIGURE  20  ftAN  FIXW  VEUXITY  VERSUS  SURFACE  NORMAL  DISTANCE  ICNDUtNSIONALIZED  BY 
(jOrPUIED  BOtfCARY  LAYER  THICMCSS  AT  CHORD  LOCATIONS  ON  UPPER  SURFACE 
(a  =  5-) 


FIGURE  21  REAN  FUX  VELOCITY  VERSUS  SURFACE  NOfiWL  DISTANCE  NONDIftNSIOttLIZED  BY 
COMPUTED  BOUCAFY  LAYER  THlCWtSS  AT  CHORD  LOCATIONS  ON  UPPER  SURFACE 
(»  =  75*) 
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FIGURE  22  ftflf)  FUW  VEIDCITY  VERSUS  SURFACE  fCPMAL  DISTANCE  NONDirtNSIOKALIZED  BY 
am/IED  BOUNDARY  WYER  THICKNESS  AT  CHORD  LjOCATICffi  ON  UPPER  SURFACE 
(o.  =  9,5-) 


FIGURE  23  l€AN  FLOW  VELOCITY  VERSUS  SURFACE  NORfW.  DISTANCE  NONDIfEMSIONALIZED  BY 
CCmnH)  BOUNDARY  LAYER  THlCWtSS  AT  CHORD  UXATIOTC  ON  UPPER  SURFACE 
(«  =  11.5') 
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APPENDIX  A 

GRID  TRANSFORyuATION'  RELATIOXSHIPS 

Several  of  the  definitions  and  coordinate  transf orr.at ion  relations 
in  the  transformeo  plane  are  summarized  below.  The  notation  used  hy 
Thames  (44)  and  Hodge  (53)  is  retained.  The  physical  plane  is  tViu  (x,  •.') 
coordinate  system  and  the  computational  plane  is  the  (  ^  ’)  coordinate 
system.  Trie  following  function  definitions  are  used; 

f(x,  y,  t)  -  a  twice  differentiable  scalar  function. 

£(>;,  y,  t)  =  i  F^(x,  y,  t)  +  j  F..,(x,  y,  t1  a  twice  d  if  f  eren  t  ia' '  ] 
vector  function  where  i  and  j  are  Cartesian  unit  vectors. 


Transforraatior.  Definitions 

^  +  y;  (A  .  1 ) 

r  =  X  y  (A.  2) 

1  =  x“-  4-  y  .  (A  .  3) 

=  (ix.,-  2r  X  +  'j>;,  ^)  (A.  4) 

P,  =  (  -..y  ,  ,  -  22  y  ,,  +  ■,>•_)  (A.  5) 

J  =  X  -V  -  X  V  ,  ( A  .  A  ) 

T  t,-  £ 

-  =  (y  -  X  p^)/J  (A.  7) 

-  =  (x^D,  -  y^;)^)/J  (A.M 

D orivatives  of  Scalar  Functions  in  Transformed  Plane 

f  =  'f/  .x  -  (v  f  .  -  V  i  )/.I  (A.-' 

f^.  =  f/.y  =  (x  f  p  (.A. 

f  =  f/.'t  =  (-f/  't)  ,  for  fixed  coordinates  (.A. ID 
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Vector  Operators  in  the  Transformed  Plane 
Gradient ; 

7f  =  '.(y^f  ,  -  y  f^)  i  +  (x  1.  -  x^  f  )  j  /J  (A.  12) 

Divergence : 

V  •  F  =  -  y  ,(F^)  +  X  iF  j  -  x^  (F,)  ,7.1  (A. 13) 

Curl : 

7  X  F  =  k-.'v  (F,,)  -  v,.(F,,)  -  X  iF  )  +  x  (F,)  7.1  (A.  14) 

-  f  -  -  ,  '  .  -1  r 

Lap  lac ian ; 

7“f  =  (.f  -  27  i  +  ('7.J")f^  +  (:/r)i  ,  (A. 15) 


Unit  Norir.al  and  langent  Vectors  in  Pi'.vsical  Plane 


N'ormal  to  '  line: 
n  -  .  /  .  ’ 


=  (-y  i  +  x-:)/' 


(A. 16) 


N'orir.al  to  f  line: 

n  ^  =  2  V  7  ■'  =  ( y 7  -  j  )  /  > 


(A. 17) 


Tangent  to  line  (increasing  v  direction!: 
-n  X  k  =  (x  jj,  +  y  _-j ) /v  'i 


(A. 18) 


Tangent  to  7  line  ^increasing  direction) : 
(  r)  f  n  "  '  '  _ 

=  -11^  X  k  =  (x^^i  +  y^j)/i  :i 


(A. 19) 


Integrals  in  the  Transformed  Plane 
Area  Integral: 

y)  d.xdy  =  S^*  f(x(t. 


.)  ,  y(  r,  r))  ,  J  d  d” 


where  R  is  the  region  R  mapped  into  the  (',  r.)  plane. 


(A. 20) 
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Vector  Line  Integral: 


F^(x,  y)  ds 


r.  ) ) .  d  r 
P 


(A. 21) 


where  contour  is  any  constant  h  line  in  the  physical  plane,  s  is  arclength 
alone  s  ,  and  h  is  the  value  of  for  the  chosen  contour  s  . 

°  ^  D 
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APPENDIX  B 


FINITE  DIFFERENCE  APPROXIMATIONS 


This  appendix  presents  the  finite  difference  approxiinat ions  that  arc- 
used  in  this  investigation.  The  differences  are  formulated  using  a  func¬ 
tion  f(-,  ’■)  in  the  computational  or  transformed  plane.  The  ’  and  ’ 
spacings  are  assumed  constant  with  value  unity.  The  truncation  error  term 
is  given  in  derivative  form  assuming  L  f  and  1’"  are  equal.  The  apprcxim.a- 
tions  are  second  order  accurate  unless  otherwise  specified.  Time  dif¬ 
ferences  are  expressed  by  Lt .  The  superscript  n  denotes  the  n^''  t  im.c 
interval  and  is  understood  when  omitted.  The  differences  are  given  for  a 
(  I,  fi)  point  location  denoted  by  subscripts  (i,  j)  and  are  understood  when 
omitted.  Space  derivatives  with  respect  to  only  are  presented  because 
the  corresponding  f  derivatives  are  identical  with  subscripts  (i,  j)  reversed. 


First  time  derivative,  backward  difference,  first  order: 
f^  =  (f"  -  f"‘^/lt  -  it  i^^/2  +  ... 

First  derivative,  backward  difference,  first  order: 

fp  =  (fj  -  fj-l)  -  fnn/2  +  ••• 

First  derivative,  central  difference: 


=  (f 


j+1 


-  f 


j-1 


)/2  -  f,_  /6  + 


First  derivative,  forward  difference: 


(B.l) 


(£.2) 


(B.3) 


(B.4) 


First  derivative,  backward  difference: 


f 


t] 


4f .  ,  +  3f .)  +  f  /3  +  . .. 
1-1  J  DDr 


(B.5) 
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Second  derivative,  central  difference: 


r" 


f  =  (f -  2f .  +  f .  -  f  /12  +  . . . 

nn  J+l  J  j-l  r:nr; 


Second  derivative,  forward  difference: 

f  =  (-f.,.  +  4f.^^  -  5f.  ,  +  2f.)  +  Ilf  /12  + 
t:',  j  +  3  j+2  j+1  2'  rr.T-, 


Second  derivative,  backward  difference: 


f  =  (f.  ^  -  4f.  ^  +  5f.  ,  -  2f.)  -  Ilf  /12  + 
Tib  1-3  1-2  1-1  2' 


Cross  derivative,  central  difference: 


^^i+1,  j  +  l  ^i+1,  j-1  ^i-1,  j-1  “  ^i-1,  i+1'^'^ 


f  ^,,)/24+  .. 


Cross  derivative,  central  in  '  and  forward  in  r  differences: 


=  -3('i+l  -  'i-l>  ^^^^+1,  i  +  1  -  ^-1,  i  +  l^ 


‘-^i+l,  i+2  “  ^i-1,  i+2^'^ 


f-,  ^ 


Cross  derivative,  central  in  f  and  backward  in  +  differences: 
^  5-  "  -^^^i+l  "  ^i-P  "  “^^^i+l,  j-1  "  ^1-1,  j-P 
^  ('1+1.  j-2  -  'i-l, 


(B.6) 


(B.7) 


(B.S) 


(B.9) 

! 


(B.IO) 


(B.ll) 
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APPENDIX  C 


GRID  MODIFICATION  FOR  WAKE  RESOLUTION 

This  appendix  describes  the  technique  that  is  used  to  insert  addi¬ 
tional  grid  lines  in  the  region  of  the  airfoil  wake.  The  body-fitted  coor¬ 
dinate  grid  system  becomes  coarse  at  distances  far  from  the  body.  The 
coarse  grid  may  reduce  the  resolution  of  the  velocity  defect  in  the  viscous 
far  wake  region.  In  order  to  increase  wake  resolution,  eight  additional  ’ 
lines  are  placed  downstream  of  the  airfoil  with  body  points  on  the  rounded 
trailing  edge.  The  original  numerically  generated  grid  has  71  '  and  UU  ^ 
lines.  Linear  interpolation  between  adjoining  7  lines  of  the  form 
=  ,x(l-l,  J)  +  x(I,  J)}/2  is  used  to  locate  a  c  line  between  lines  1=1 
and  I  =  2,  I  =  2  and  1  =  3 ,  I  =  69  and  1  =  70,  and  I  =  70  and  I  =  71.  The 
procedure  is  repeated  using  the  new  lines  as  well  except  only  one  additional 
line  is  located  below  the  1=1  cut.  In  this  way  the  fine  grid  spacing  has 
an  asymmetry  for  better  wake  resolution  at  angle  of  attack.  The  final  79 
by  44  grid  is  shown  in  Figure  3.  Table  C.I  gives  the  numbered  1  index 
(F  line)  designation  for  the  grid  systems  in  the  wake  region. 
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APPENDIX  D 


CONTOUR  INTEGRAL  METHOD  FOR  DETERMINING  BODY  FORCE  COEFFICIENTS 

The  determination  of  the  resultant  force  that  a  flow  field  exerts  on 

a  body  usually  involves  the  calculation  of  surface  pressure  and  viscous 

stresses  and  the  summation  of  these  forces  over  the  body  surface.  An 
alternate  approach  is  to  apply  a  control  volume  analysis  to  a  region 
enclosing  the  body.  In  the  foregoing  discussion,  a  control  volume  analysis 
for  a  body  in  a  two  dimensional  flow  is  presented. 

Consider  a  general  body  immersed  in  a  flow  field  with  surface 
Define  a  fixed  control  volume  V  with  outer  surface  and  inner  surface  C-^. 

Conservation  of  linear  momentum  for  the  fluid  in  V  is  expressed  by  the 

Cauchy  Integral  Equation  of  Motion 


(D.l) 


where  .  is  the  density,  v  velocity  field,  ^  body  force  per  unit  mass,  n 
unit  outwara  surface  normal  vector,  T  traction  stress,  and  C  is  the  total 
surface  Apply  the  Reynolds  Transport  Theorem  to  the  material 

derivative  term  and  the  Divergence  Theorem  to  the  surface  integral  and  obtain 
from  Equation  D.l  for  the  i  Cartesian  vector  component 

/  V —  (,  V  .  )  dV  +  I  c  ,  (.  V  .  V .  -  T  .  .  )  dV  =  0  ( D .  2 ) 

V  \' 

Applv  the  Divergence  Theorem  to  the  second  volume  integral  and  express  tin- 
resulting  surface  integral  as  separate  integrals  over  and  \  to  obtain 
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/4—  (p  V  . )  dV  +  I  n  .  (p'  V .  V .  -  T  .  . )  dc  +  /  n  .  (.  v  .  v  .  -  T  .  . )  d  ' 

"t  1  vy  J  1  J  Ji  J  1  J  Ji 


(D.3) 


Now,  on  the  solid  non-porous  bodv  surface  n.  v.  =  0.  The  integral 

B  j  j 

of  the  resultant  traction  force  over  c  is  the  net  force  exerted  bv  the 

D 

body  on  the  control  volume  fluid.  Therefore,  the  force  exerted  by  the 
fluid  on  the  body  la  given  by 


/  "j 


T.  .  dc 


Substitute  Equation  D.4  into  D.3  and  use  the  above  surface  condition  to 
obtain  a  general  expression  for  the  force  exerted  by  the  fluid  on  the  body. 


' .  =  -  I  r  V  .  dV  +  /  n  .  (T  .  .  -  V  .  V  , ' 

1  ^t  J  1  Jjji  ij 


Now  consider  the  case  of  two-dimensional  flow  in  the  x-y  plane.  Then 
Equation  D.5  becomes  the  following  expression  for  the  force  per  length  of 


span  ^ 


.  =  -  T--  /  .  V  .  dA  +  /  n ,  (T .  .  -  .  V  . 

1  t  J  1  J  J  Ji  1 


V  , )  ds 
J 


where  A  is  the  cross  sectional  area  of  the  control  volume  and  s,  is  the  outt 

perimeter  of  A  as  shown  in  Figure  36. 

Substitute  for  the  traction  stress  1..  =  -p'\  .  +  in  Knuaticni  I’.n 

I  1  J 1  J 1 


and  obtain 


=  /  n  .  (-p -  v.v.ids  -  7 —  /  .  V .  dA 

^  J  J  Ji  Ji  ij  tj  1 


where  p  is  the  pressure,  ' is  the  laminar  viscous  stress  auJ 
is  the  Kronecker  Delta  Function.  Next  express  the  (low  variahles  in 
terms  of  turbulent  mean  and  fluctuating  variibies,  assume  incompre.- s  i ;  1 1 
flow,  and  Reynolds  average  Equation  0.7  to  obtain 


c  =  /  n. 

>■  .  J 


-  p  .  .  +  .  .  -  (  V  .  V  .  +  V  !  V  1  )  d  s 
Ji  jl  ^  i  1  J 


-  -t  f  \ 


The  term  -  v!  v!  is  the  turbulent  Kevnolds  stress  ...  let  thi,  tele’,  stre 
1-  J  ■  i  J  1 

be  defined  as  ' .  .  =  '  .  .  +  \  .  ..  Then  Equation  D .  E  becor.es 
I  J  1  1 1  f  J  1 


1  ^  1 


/  n  .  -p  ■  ,  .  +  .  v  .  V  .  us 

'  J  J  r  t  ,t  1  1  .1 


Introduce  the  following  nond  imensional  variables  in  roueition 


P  -  P.. 


\' . 
1 

Vi  = 


and  t  =  —  ,  and  in  tt  id,.; 


Iji  ■  V,.  (Ji 


where  ;■  is  the  fieestiear.  VL-locitv  and  !  is  a  riliis 


length  in  tlie  x  direction.  Also,  define  thi'  forci-  eoi’f  t  ic  i  ent  =  'i 

i  1  ' 


and  obtain  from  Equation  1'.'^ 
(\  =2  S  n  .  -p  •  ,  ,  +  , 


^  =2  J  n  .  -  p  ■  .  .  +  ,T  ,  V .  .  d  s  - 

f  .  J  '  ,!■  Iji  >  J 


where  all  variables  art  understood  to  be  mean  d  im.t  ns  io:i  1  ess  v.iri; 

the  turbulent  eddy  viscosity  concept  is  used  whi-re 

.  .  =  (.  +  •  ,1  •  (  •  .V,  +  •  ,v.  ) 

I  J  I  >•  1  J  J  1 


and  define  He^  =  ,  - 


which  transforms  Ecuation  D .  10  te  the  'el  lav,  in, 


0.  =2 


-p  '  .  ,  +  -  (  . V .  +  . V  )  -  V . V .  ds  -  2  - 


(0.11  I 


Expand  Equation  D.ll  into  x  and  y  compcinent  equations 


I 


1 

Re. 


2.,  ^  + 
1  x 


(  -  +  1  -  (u  n, 

.•V  -x  1 


+  u\’n„)  ds 


C  . 

i 

V 


(D.12a) 


(uvr.j  +  Vn ,, )  ds 


(D.12b) 


w’acru  iij  ar.d  n,,  ari-  Ihi.-  x  and  y  components  of  the  unit  outward  norrail 
Vector  n  for  t!ie  outer  bovindary  s_  of  area  A.  I'quations  (!?.12a'>  and  ' 

are  expressions  for  tJie  x  and  y  force  coefficients  exerted  by  a  two- 
dimensional  Inconpressihle  turbulent  flow  on  an  arbitrary  two-d  ir-.t/ns  i  ona  1 
body  wliere  path  s,  encloses  the  body  and  an  eddy  v.scosity  irathod  is  used, 
to  model  turbulence. 


Lift  and  draj:  coefficients  are  then  obtained  from  the-  force  cc-ef  1  ic  it-nt 
as  folK'ws  wheri-  >  is  the  geometric  angle  of  atlacL. 


ij  =  Cj,  cos  I  -  Cj.  sin 

‘  V  ‘  X 


( b .  1  n 


L,  =  (  .  Sin  :  +  r  cOi 

n  ;  t 

\’  X 


I  L.  h.  1 


i1 


APPKNI'IX  F. 


TI^EJNCRFX-l^Er.T  INFLI'ENCF  ON  SOR  CONA^^FROFFXF 

Implicit  SUR  methods  have  been  shuvn  by  linear  stability  analytes 
to  be  unrond  i  t  iona  1 1  y  stable.  lU'Wevt  r  ,  convergence  has  required  a  sr.a  1  I 
time  step  size  in  tliis  investigation  and  elsewhere.  The  following  analy¬ 
sis  indicates  the  relationship  between  the  time  step  size  and  cemvergenc e . 

Consider  the  systc:"  of  difference  eqiiatiens  to  be  writteri  i:;  tin 
following  matrix  iort;: 


.•\v  =  b  ( i: .  1 ) 

where  .A  is  the  coefficient  matrix,  F  is  thr.-  vector  of  unknown  prit.it  ivi- 
variables  at  time  step  n  with  u,  v,  p  grouped  together  for  each  (i,;l 
point  location,  and  b  is  a  vector  oi  "constants".  For  example,  the  x  cor 
ponent  of  the  momentum  equation,  Fquatiem  92,  at  the  (i.,j)  point  is 

Q  u^^^  =  o'? .  ^  +  u‘.t[-YKTA  (p.,,  .  -  p.  ,  .)  +  YXI  (p,  -  P.  .  1 1 

ij  ij  1  +  ]  J  1-1  .)  1  J  +  1  1  J-l 


+  rc  (.1.U 

ICl 

.  -  u 

J 

IC2 

^v:  ^ 

(4u. 

,TV1 

+  c\yj\  (u. 

j+i 

+  u  . 
1 

j-1 

1'  ■  1  •  n  Y  ' 

1-1  j+i 

-  u 

(s-1) 

ij 

n 

,)  +  VC  (4u, 


i  dCl 
1 


“i  JV2^  RFT  ^"i+1  j  ^'i-1  ,i  ^ 


i+1  i+1 


^i+1  i-1  ‘‘i-1  i-1 


d'.:) 


where  all  coefficients  are  evaluated  at  point  (i,j)  and  prit.itivc  variab 
are  at  time  step  n  unless  given  otherwise,  and  n  is  given  by 

Q  =  1  +  3  .‘.tt  rc 

Rearrange  Equation  F..1  and  obtain 


VC  +  UV  +  V\  )  +  - —  (  ALFA  +  ) 

Kf 

I 


'C  u 


ic:  j 


.'V  ,  Cl 


■-  J 


^  BETA 

^  reY  "i- 


.  ,  +  4  IT  n  ,  ,  .  +  4  'dv 

.1-1  icl  J 


IVl  i 


t  4i|i  „  ,  +  YEiA  ..  ,  ,,,  -  VC  u.  -  VV 

RET  1-1  j  1-1  J  RET  1-1  J+1  1  JCJ  1  J\. 


+  u.  .  +4  VC  u.  ,  +4  VV .  u  .  ,  -  YXI  p.  .  .  -  ^ 

RET  1  j-1  '  ‘  1  JCl  '  ■  1  JVl  ‘  1  j-1  Ll  ij 

,  CAMA  BETA  ^ALFA  . 

RKl  j  +  1  j+1  ~  RET  “i+1  j-1  RET  “i+1  j  “  '’i+1  j 


BETA 

RET  ’  ^‘i+1  j  +  1 


u*?  .  ^  (s-1) 

-  _i^  +  u  D 


A  sir.ilar  y  component  of  rnomentur'.  equation  is  obtained  if  u  is  rep-laced  ly 
v;  -  YETA  by  xrrA;  and  YXI  by  -  XX] . 

T!je  finite  difference  pressure  Equation  94  becones 

-  TV  p,,..,  .  +  BETA  p.  .  ,  +4  TV  p  ,  .  +  ALFA  p.  ,  . 

^r,  j  J  1-1  J-1  ‘  I\  1  J  '  1-1  J 

-  BETA  P.  ,  -  VV  p.  +  G.A'AA  p.  .  ,+4  VV  p,  -Bp,. 

1-1  J  +  1  ‘  1  J\2  ‘  1  J-1  '  1  .1\  1  '  ij 

+  G.V14  p.  -  BETA  p._^,  .  ,  +  ALFA  p._,,  .  +  BETA  p._,,  ,  =  -K:if 

1  j  +  1  1+1  J-1  '^1+1  J  '  1  +  ]  j+  1 


where  B=2  (ALFA+G/ALA)  +  1  (  rv,+  VV  )  and  RHS  =  -  +  [  (F.X)  "+2  (VX)  (I’Y )+ iVY  » ‘  ' 

Convergence  ’'heorcms  (SI)  for  systems  of  linear  equations  wiili  con.^tant 
coefficients  show  that  a  necessary  and  sufficient  condition  for  convergence 
requires  that  the  coefficient  matrix  be  diagonally  dominant.  If  tb.i.--  criteria 
is  applied  to  the  above  quasi-linear  system,  the  following  convergence  criteria, 
inequality  E.b  from  the  momentum  equation  and  inerjualitv  E.y  f  i\im  t’le  pt't.^se.ri. 


equation,  are  obtained: 

5  (  rc  +  VC  +  rv  +  vv  )  +  (  at  fa 

+  BETA  +  2(  YXI  +  YEIA  1  •  9 

K  h  1  ■  t 


+  (;a:t\  ' 


or  ■  J  (  1 'C  +  VC 


which  becomes 


rv  +  VV  )  +  ;  (  ,.\lfa  +  ga-v,  •  + 

K  r.  1  f 


.'.t-  2  (  IT  +  VC  +  rv  +  vv  )  +  beta  +.'1 


f !  .  '■  I 


and 


2  (  i  ALFA  ,  +  CAM.\  '  )  +A  ,  BETA  +  5  (  ,  I'V  +  VV  ,  ')  • 

2(;aLFA,  +  ,GA.'L-\  )  +  3  (;L’V  +  '  VV ,  )  (E.7) 

or  A  ;bETa'  +  2  ( ' TV ' +  VV  )  ^  0  (E.8) 

The  convergence  criteria,  Equation  E.6,  provides  a  lirr.it  on  It  a 
function  of  both  the  grid  geometry  and  flow  field  solutio;;.  Ihu  rig'nt  hand 
side  of  the  inequality  has  its  smallest  magnitude  near  tb.e  trailing  edge 
and  increases  rapidly  away  from  the  airfoil.  An  order  of  magnitude  analysis 
of  the  grid  coefficients  and  an  examination  of  the  computational  solutions 
indicate  that  diagonal  dominance  occurs  if  '.t  <  0.0003  for  points  near  tie 
trailing  edge,  L  '  0.001  elsewhere  near  the  body,  and  t  ■  1  in  tlu-  far 

field.  The  pressure  equation  convergence  criteria.  Equation  E . £ ,  has  no  t 
dependence  and  Is  approximately  satisfied  only  at  points  far  from,  the 
body  where  the  left  hand  side  approaches  zero.  These  criteria  can  only 
indicate  a  possible  time  step  restriction  since  the  system  of  equations  l;as 
variable  coefficients  and  lower  order  nonlinear  it ies . 

A  numerical  experiment  was  conducted  to  determine  the  time  step  cfftct 

on  convergence.  A  computation  for  the  flow  over  the  NACA  0012  airfoil  at 

an  angle  of  attack  of  eight  degrees  was  used.  The  solution  was  advanced 

one  time  step  with  the  time  step  size  and  number  of  iterations  per  time 

step  varied.  The  solution  converged  for  time  steps  of  0.0005  and  0.001  while 

diverging  for  a  time  step  of  0.005.  The  relative  maximum  error  magnitudes 

for  the  u  component  of  velocity,  surface  pressure  ,  and  field  pre.ssure 

p  as  functions  of  number  of  iterations  are  shown  for  each  t imt  step  in 

Figures  37,  38,  and  39.  The  relative  error  is  defined  as  (f  -  f  .''/f  . 

n  n-1  n 

The  V  component  of  velocity  errors  behave  .similarlv  to  tb.t  u  cor.penent. 

Next,  the  pre.s.suri  was  held  fixed  and  the  convergence  of  the  \-el  oc  it  v 
was  monitored  (Figures  37  ,  38,  and  I'Vi  .  qhe  rapid  rate  e'f  ronvt  rgeiici  i 


observed.  The  u  and  v  velocity  components  were  also  held  fixed  i 
turn  to  test  convergence.  In  both  cases,  convergence  was  similar 
the  general  solution  case,  These  numerical  experiments  indicate  t 
the  convergence  of  pressure  is  slow  and  that  convergence  is  depen 
on  time  step  sizes  similar  to  those  predicted  by  the  simple  analy 
abo\-e  . 


APPE;;rjix  i- 

FAR  FIELD  BOU.VDARY  CONDITIONS 

TIk-  far  field  boundary  for  a  nunerical  solution  must  usually 
bf  at  a  finite  distance  from  the  body  of  interest.  This  constraint 
may  pose  difficulties  if  the  only  known  far  field  boundary  conditions 
are  those  at  infinity.  This  analysis  uses  complex  variable  methods 
for  incompressible  potential  flow  and  develops  approximate  far  field 
boundary  conditions  for  use  at  a  finite  distance. 

Consider  a  circular  cylinder  radius  a  and  center  at  the  cricin. ,  wit':, 
positive  clockwise  circulation  I,  in  a  uniform  stream  (velocity  V  1 , 
at  angle  of  attack  This  coordinate  system  is  identical  with  the 
physical  plane  of  the  airfoil.  The  complex  velocity  is 

)  i.t  _ 

w  (z)  =  r..  e  -  — —  -7—:: —  +  i  —  “  (F .  11 

Z  “  z 

where  the  first  term  is  the  uniform  stream,  contribution,  the  second 
term  is  a  doublet  of  strength  2”a'U  ,  and  the  third  term,  is  a  vortex  of 
strength  7. 

N'e.xt ,  use  the  Kutta-.Joukowski  Theorem  ^or  lift  and  the  definition 
of  the  lift  coefficient  to  ob  ta  in 

:■  =  (1/2)  Cj^l'^c  (F.2) 

where  c  is  the  airfoil  chord. 

Nondimensional  ize  Equation  F.l  velocities  with  the  freestream  velocit>- 
and  lengths  witli  the  airfoil  chord  c  and  substitute  Equation  1.2  into  F.l 


and  obtain 


w  (Z)  =  -  (2)' 


+  i 


I.  1 


(F.  3) 


where  U’  is  the  nondimensional  complex  velocity  and  ?.  is  the  nondimensional 
complex  variable  Z  =  x  +  iy. 

Airfoil  transformations  such  as  the  Joukowski  and  Treffetz  Trans¬ 
formations  transform  a  cylinder  into  an  airfoil  witli  a  chord  of  four 
radii.  Thus,  for  a  given  chord  c,  choose  the  cylinder  radius  to  he  1/4 
c  which  gives  a  scaling  factor  ^  =  1/4  in  the  doublet.  Then  the  complex 
Velocity  for  the  flow  around  the  scaled  cylinder  becomes 


C 


W  (Z)  =  e" 


+  1 


1.  1 


(F.4) 


16Z'  '  '4  Z 

Expression  F.4  can  then  be  used  to  approximate  the  far  field  potential 
flow  over  the  airfoil  of  chord  c  with  circulation  T  in  the  same  physical 
plane . 

Let  the  far  field  boundary  be  a  circle  of  radius  where 


Z  =  r^  e  - ,  0  ■_ 


Substitute  for  Z  in  Equation  F.4  and  obtain 
C, 


kj.  (  )  =  e 


Ihr  .■ 
t 


+  i^  -  e- 
"f 


(F.5) 


Combine  the  exponential  terms  in  F.5  and  use  the  fact  that  V  =  u-iv, 

where  u  and  v  are  the  x  and  y  velocity  components,  to  obtain 

C, 


/  N  -  1 

U,(  )  =  COS  (  -  — 

f  16r . 


1 


cos  (.1  -  2C)  +  ycr  sin 
4  r  - 


t 


C. 


V'  (  )  =  sin  a  +  i  sin  (.<  -  2:)  -  7—  cos 


(F.bel 

(F.6b) 


In  both  Equations  F.6,  the  first  term  is  the  freestream  velocity  contribution, 
the  second  term  is  the  doublet  contribution,  and  the  third  term  is  tht' 
vortex  contribution.  As  r^.  approaches  infinity,  the  infinity  freestream 
boundary  conditions  for  u  and  v  are  recovered.  Alsti,  tlie  doublet  terr: 
is  of  higher  order  in  1/r^  than  the  vortex  term  for  the  case  of  n  lifting 


airfoil  and  mav  be 


d. 


1 10 
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The  far  field  pressure  is  obtained  by  applying  Bernoulli's 


equation  for  irrotational  incompressible  flow  and  using  the  nopdimensional 
fo 'm  for  pressure.  The  result  for  the  far  field  boundary  is 

Pj  C)  =  1/2  {l-Cu^-  4.  Vj,-')'-  (F.7) 

The  far  field  boundary  conditions  of  Equations  F.6  and  r.7  modify  th.e 
infinity  freestream  conditions  by  incorporat ing  effects  on  the  flow  of 
body  thickness  and  circulation  at  a  large  but  finite  distance  from 
the  airfoil. 

Transonic  small  disturbance  theory  for  slender  bodies  and  airfoils 
has  similar  expressions  for  the  far  field  conditions.  Klunker  (101)  has 
obtained  an  asymptotic  solution  applicable  at  large  distances  from,  thin 
airfoils.  This  form  has  been  used  as  a  numerical  outer  boundary  condition. 
For  the  limiting  case  of  Incompressible  two-dimensional  flow,  the  doublet 
and  vortex  strengths  are  compared  below  with  the  general  form.s  found  in 
Equation  F.l. 

The  nondimensional  velocity  potential  doublet  in  the  far  field 
from  Klunker  becomes 


•d  '  2-  2/F(s)ds  (F.8) 

where  F(s)  is  the  airfoil  thickness  function.  The  doublet  sfength  is  thu.- 

2 ^F(s)ds.  The  NACA  four  digit  airfoil  thickness  function  (A7)  is 

F(s)  =  (.2969  (-)*'  -  .126  -  -  .3516(-)  +  .2843  (-)  - 

0.2  c  c  c  c 

.1015  (-)  )  (F.Q) 

c 

where  c  is  the  airfoil  chord,  f  is  the  maximum  thickness  as  a  fraction 
of  chord  and  s  is  the  chord  distance,  0  ^  s  ^  c. 

Substitute  the  thickness  function  Equation  F.9  into  the  doublet  strengtli 
expression  and  obtain  the  following  doublet  strength  from  transonic  sm.all 


disturbance  theorv; 


0.685  fC 


(!■.  10) 
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For  the  case  of  a  NACA  0012  airfoil,  hp  =  0.0822c‘.  The  corresponding 

doublet  strength  in  the  simple  potential  flow  model  from  Equation  F.l  is  T' 

c  ^  ^ 

(^)  or  .3972c  .  Klunker  points  out  that  the  doublet  may  be  neglected  in 

the  far  field  for  finite  circulation  flows  because  the  vortex  contribution 

1  1  ^ 

is  of  order  (—  )  compared  with  an  order  (—  )  for  the  doublet. 

’^f  ^^f 

The  vortex  potential  term  of  Klunker  for  the  case  of  incompressible 
two-dimensional  flow  is 

Vy  sgn  (y)  +  tan  ^  (^) '■  (F.ll) 

where  F  is  the  clockwise  circulation,  the  inverse  tangent  function  is 
defined  on  the  interval  (-  ^  ■^)  ,  and  the  angle  within  braces  is  defined 
clockwise  from  the  negative  x  axis.  This  potential  function  is  identical 
with  the  vortex  function  in  Equation  F.l  where  the  clockwise  strength  s  again 
r  and  the  angular  orientation  is  the  right  hand  polar  coordinate  dei  ed 
counterclockwise  from  the  positive  x  axis. 

Thus  the  far  field  potential  for  each  method  differs  only  in  the 
thickness  or  doublet  term.  Furthermore,  for  flows  with  finite  cir'-  ration 
this  difference  is  negligible  because  the  vortex  term  dominates. 
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APPENDIX  G 


NUMERICAL  INVISCID  FLOW  RESULTS 


A  numerical  solution  for  the  inviscid  flow  over  the  NACA  0012 
airfoil  is  obtained  for  comparison  with  both  the  Kaiver-Stokes  computation 
and  the  experimental  data.  The  numerical  algorithm,  developed  by  Thames 
(44),  which  solves  Laplace's  equation  for  the  stream  function  in  two- 
dimensional  flow  is  used.  The  method  uses  SOR  iteration  to  solve  the 
system  of  equations  written  in  terms  of  body-fitted  coordinates.  Tlie 
stream  function  value  on  the  body  can  either  be  specified  or  computed  by 
imposing  a  given  circulation  or  a  Kutta  condition. 

The  stream  function  equation  for  two-dimensional  flow  is 


V  +  (  =0  (G.l) 

XX  yy 

Rewrite  Equation  G.l  in  terms  of  the  transformed  variables  in  the 
computational  plane  using  expressions  in  Appendix  A  to  obtain 


y  r  i 


+a  li,- 


+  av,.  +T  4'  j  =  0 


((■..2) 


Fr,  ■■  "nri 

The  boundary  condition  for  the  stream  function  on  the  outer  boundar;.  ,  wb  L' 
is  a  radius  of  ten  chords,  is  the  freestrcam  value 


y^  cos  ij  -  Xj  sin 


(1  .  M  ' 


where  (x^,  y^)  are  coordinates  on  the  outer  computational  boundary  a  ' 
is  the  angle  of  attack.  Hodge  and  Ghia  (82)  have  shown  tiy  an 
inviscid  analysis  that  this  boundary  condition  induces  an  error  in  lift 
of  less  than  1%  for  angles  of  attack  less  than  iO'’.  The  boundary 
condition  on  the  body  is 

V,  =  constant  (('..3b' 

b 

The  constant  can  be  specified  or  determined  by  imposing  a  given  circula¬ 
tion  or  a  Kutta  condition.  The  Kutta  condition  which  matches  the  upper 
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and  lower  trailing  edge  surface  velocities  by  extrapolation  is  select¬ 


ed.  The  details  for  each  option  are  given  by  Thar'.es  (44). 

The  derivatives  in  the  transfcrmed  Equation  0.2  are  approxir.at ed 
by  central  difference  formulas  found  in  Appendix  B.  The  resulting 
system  of  equations  is  solved  by  SOR  as  described  by  Tl.ame^  (44  i.  S  lo¬ 
tions  are  obtained  for  0°,  5°,  7.5°,  9.5°,  and  11.5'‘  for  the  h'AO.'i  14 
airfoil  section.  The  body  surface  pressure  coefficients  can  be  calculated 
from  the  stream  function  solution  as  follows.  The  body  pressure  coefficie 
is  defined  as 

C  =  1  -  (u-  +  V")  (i  .4) 

P 

where  at  a  location  on  the  body  u  and  v  are  the  x  and  y  component.s  of  the 
flow  velocity  nondimensional  ized  by  the  freestrear.  velocity.  Ihse-  the 
definition  of  the  stream  function  to  obtain 


Then  rewrite  the  Equation  C.5  in  terms  of  the  transformed  variables  cind 
note  that  1  =  0  on  the  body  to  give 

C  =  1  -  %  (■:  )-  (G.b) 

p  J  n 

where  all  quantities  are  evaluated  on  the  body  surface.  The  body  pressure 
coefficients,  calculated  using  Equation  G.6  with  second  order  one-sided 
differences,  are  plotted  in  Figures  16,  17,  18,  19,  and  33.  A  typical 
streamline  plot  for  A  =  11.5°  is  shown  in  Figure  AO.  The  numerical  in- 
viscid  solution  is  also  compared  with  an  analytical  solution  (67)  which 
uses  Theodorsen's  (12)  method.  The  excellent  agreement  is  seen  in 
Figure  AO  where  both  solutions  are  plotted. 
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APPENDIX  H 


EXPERIMENTAL  DATA 


This  appendix  presents  the  unpublished  data  which  have  been  obtained 
experimentally  at  the  Air  Force  Frank  J.  Seiler  Research  Laboratory,  FSAF 
Academy,  Colorado,  by  laboratory  personnel  and  used  in  this  work. 

Airfoil  upper  and  lower  surface  mean  pressure  coefficient  measure¬ 
ments  for  nominal  attack  angles  of  0°,  5“,  8°,  10°,  and  12°  at  a  chord 
Reynolds  number  of  170,000  are  given  in  Table  H.l.  The  pressure  coefficient 
is  defined  as  the  static  pressure,  relative  to  freestream  static  pressure, 
nondlmensional ized  bv  tlie  freestream  dynamic  pressure  as  follows:  C  = 


'p  -  P,)/.  ■  l',7  •  sets  of  hot  wire  anenometry  velocity  field  data 

near  the  upper  and  lower  airfoil  surfaces  on  designated  paths  are  presented 
in  Table  h.ll.  The  velocity  field  and  Reynolds  stress  u'v'  data  for  tlie 
near  wake  region  arc  presented  in  Table  H.III.  The  data  have  been  non- 
dimensional  i  7,ed  by  the  measured  freestream  velocity.  The  spatial  locations 
in  the  physical  coordinate  system  for  each  path  are  given  in  Table  H . IV 
where  location  (D,0)  is  the  airfoil  center  and  the  x  axis  lies  along  the 
airfoil  chord. 

The  spatial  measurement  locations  correspond  with  computational  grid 
point  locations  on  constant  ^  lines  or  '  lines  in  the  physical  plane.  The 
experimental  lead  time  necessitated  the  use  of  the  original  71  >;  44  grid 
system.  The  I,  J  index  notation  designates,  however,  the  point  indexing 
(•','■)  in  the  final  79  x  44  grid  described  in  Appendix  C. 
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TABLE  [[.I  Experimental  Surface  Pressure  Coefficients  for  NACA  0012 
At  Various  Angles  of  Attack,  Re  =  170,000 
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TAKI.K  11. IT  Kxpcrimontal  Wlocity  Field  Measurements  (Two  Sels)  Near  SurFacc  of  NACA  0012  Airfoil 
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